Home

Context Navigation

source: palm/trunk/SOURCE/flow_statistics.f90 @ 1652

Last change on this file since 1652 was 1652, checked in by raasch, 9 years ago
bugfix in calculation of vertical profile of energy production by turbulent transport of TKE bugfix for output of configuration settings (mbuild) output of login init command (mbuild)
Property svn:keywords set to `Id`
File size: 143.0 KB

Line
1	#if ! defined( __openacc )
2	SUBROUTINE flow_statistics
3
4	!--------------------------------------------------------------------------------!
5	! This file is part of PALM.
6	!
7	! PALM is free software: you can redistribute it and/or modify it under the terms
8	! of the GNU General Public License as published by the Free Software Foundation,
9	! either version 3 of the License, or (at your option) any later version.
10	!
11	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
12	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
13	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
14	!
15	! You should have received a copy of the GNU General Public License along with
16	! PALM. If not, see <http://www.gnu.org/licenses/>.
17	!
18	! Copyright 1997-2014 Leibniz Universitaet Hannover
19	!--------------------------------------------------------------------------------!
20	!
21	! Current revisions:
22	! -----------------
23	! bugfix in calculation of energy production by turbulent transport of TKE
24	!
25	! Former revisions:
26	! -----------------
27	! $Id: flow_statistics.f90 1652 2015-09-17 08:12:24Z raasch $
28	!
29	! 1593 2015-05-16 13:58:02Z raasch
30	! FORTRAN errors removed from openacc branch
31	!
32	! 1585 2015-04-30 07:05:52Z maronga
33	! Added output of timeseries and profiles for RRTMG
34	!
35	! 1571 2015-03-12 16:12:49Z maronga
36	! Bugfix: output of rad_net and rad_sw_in
37	!
38	! 1567 2015-03-10 17:57:55Z suehring
39	! Reverse modifications made for monotonic limiter.
40	!
41	! 1557 2015-03-05 16:43:04Z suehring
42	! Adjustments for monotonic limiter
43	!
44	! 1555 2015-03-04 17:44:27Z maronga
45	! Added output of r_a and r_s.
46	!
47	! 1551 2015-03-03 14:18:16Z maronga
48	! Added suppport for land surface model and radiation model output.
49	!
50	! 1498 2014-12-03 14:09:51Z suehring
51	! Comments added
52	!
53	! 1482 2014-10-18 12:34:45Z raasch
54	! missing ngp_sums_ls added in accelerator version
55	!
56	! 1450 2014-08-21 07:31:51Z heinze
57	! bugfix: calculate fac only for simulated_time >= 0.0
58	!
59	! 1396 2014-05-06 13:37:41Z raasch
60	! bugfix: "copyin" replaced by "update device" in openacc-branch
61	!
62	! 1386 2014-05-05 13:55:30Z boeske
63	! bugfix: simulated time before the last timestep is needed for the correct
64	! calculation of the profiles of large scale forcing tendencies
65	!
66	! 1382 2014-04-30 12:15:41Z boeske
67	! Renamed variables which store large scale forcing tendencies
68	! pt_lsa -> td_lsa_lpt, pt_subs -> td_sub_lpt,
69	! q_lsa -> td_lsa_q, q_subs -> td_sub_q,
70	! added Neumann boundary conditions for profile data output of large scale
71	! advection and subsidence terms at nzt+1
72	!
73	! 1374 2014-04-25 12:55:07Z raasch
74	! bugfix: syntax errors removed from openacc-branch
75	! missing variables added to ONLY-lists
76	!
77	! 1365 2014-04-22 15:03:56Z boeske
78	! Output of large scale advection, large scale subsidence and nudging tendencies
79	! +sums_ls_l, ngp_sums_ls, use_subsidence_tendencies
80	!
81	! 1353 2014-04-08 15:21:23Z heinze
82	! REAL constants provided with KIND-attribute
83	!
84	! 1322 2014-03-20 16:38:49Z raasch
85	! REAL constants defined as wp-kind
86	!
87	! 1320 2014-03-20 08:40:49Z raasch
88	! ONLY-attribute added to USE-statements,
89	! kind-parameters added to all INTEGER and REAL declaration statements,
90	! kinds are defined in new module kinds,
91	! revision history before 2012 removed,
92	! comment fields (!:) to be used for variable explanations added to
93	! all variable declaration statements
94	!
95	! 1299 2014-03-06 13:15:21Z heinze
96	! Output of large scale vertical velocity w_subs
97	!
98	! 1257 2013-11-08 15:18:40Z raasch
99	! openacc "end parallel" replaced by "end parallel loop"
100	!
101	! 1241 2013-10-30 11:36:58Z heinze
102	! Output of ug and vg
103	!
104	! 1221 2013-09-10 08:59:13Z raasch
105	! ported for openACC in separate #else branch
106	!
107	! 1179 2013-06-14 05:57:58Z raasch
108	! comment for profile 77 added
109	!
110	! 1115 2013-03-26 18:16:16Z hoffmann
111	! ql is calculated by calc_liquid_water_content
112	!
113	! 1111 2013-03-08 23:54:10Z raasch
114	! openACC directive added
115	!
116	! 1053 2012-11-13 17:11:03Z hoffmann
117	! additions for two-moment cloud physics scheme:
118	! +nr, qr, qc, prr
119	!
120	! 1036 2012-10-22 13:43:42Z raasch
121	! code put under GPL (PALM 3.9)
122	!
123	! 1007 2012-09-19 14:30:36Z franke
124	! Calculation of buoyancy flux for humidity in case of WS-scheme is now using
125	! turbulent fluxes of WS-scheme
126	! Bugfix: Calculation of subgridscale buoyancy flux for humidity and cloud
127	! droplets at nzb and nzt added
128	!
129	! 801 2012-01-10 17:30:36Z suehring
130	! Calculation of turbulent fluxes in advec_ws is now thread-safe.
131	!
132	! Revision 1.1 1997/08/11 06:15:17 raasch
133	! Initial revision
134	!
135	!
136	! Description:
137	! ------------
138	! Compute average profiles and further average flow quantities for the different
139	! user-defined (sub-)regions. The region indexed 0 is the total model domain.
140	!
141	! NOTE: For simplicity, nzb_s_inner and nzb_diff_s_inner are being used as a
142	! ---- lower vertical index for k-loops for all variables, although strictly
143	! speaking the k-loops would have to be split up according to the staggered
144	! grid. However, this implies no error since staggered velocity components are
145	! zero at the walls and inside buildings.
146	!------------------------------------------------------------------------------!
147
148	USE arrays_3d, &
149	ONLY: ddzu, ddzw, e, hyp, km, kh, nr, p, prho, pt, q, qc, ql, qr, qs, &
150	qsws, qswst, rho, sa, saswsb, saswst, shf, td_lsa_lpt, td_lsa_q,&
151	td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us, usws, &
152	uswst, vsws, v, vg, vpt, vswst, w, w_subs, zw
153
154	USE cloud_parameters, &
155	ONLY: l_d_cp, prr, pt_d_t
156
157	USE control_parameters, &
158	ONLY: average_count_pr, cloud_droplets, cloud_physics, do_sum, &
159	dt_3d, g, humidity, icloud_scheme, kappa, large_scale_forcing, &
160	large_scale_subsidence, max_pr_user, message_string, ocean, &
161	passive_scalar, precipitation, simulated_time, &
162	use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
163	ws_scheme_mom, ws_scheme_sca
164
165	USE cpulog, &
166	ONLY: cpu_log, log_point
167
168	USE grid_variables, &
169	ONLY: ddx, ddy
170
171	USE indices, &
172	ONLY: ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums, &
173	ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzb_diff_s_inner, &
174	nzb_s_inner, nzt, nzt_diff
175
176	USE kinds
177
178	USE land_surface_model_mod, &
179	ONLY: dots_soil, ghf_eb, land_surface, m_soil, nzb_soil, nzt_soil, &
180	qsws_eb, qsws_liq_eb, qsws_soil_eb, qsws_veg_eb, r_a, r_s, &
181	shf_eb, t_soil
182
183	USE pegrid
184
185	USE radiation_model_mod, &
186	ONLY: dots_rad, radiation, radiation_scheme, rad_net, &
187	rad_lw_in, rad_lw_out, rad_sw_in, rad_sw_out
188
189	#if defined ( __rrtmg )
190	USE radiation_model_mod, &
191	ONLY: rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
192	#endif
193
194	USE statistics
195
196	IMPLICIT NONE
197
198	INTEGER(iwp) :: i !:
199	INTEGER(iwp) :: j !:
200	INTEGER(iwp) :: k !:
201	INTEGER(iwp) :: nt !:
202	INTEGER(iwp) :: omp_get_thread_num !:
203	INTEGER(iwp) :: sr !:
204	INTEGER(iwp) :: tn !:
205
206	LOGICAL :: first !:
207
208	REAL(wp) :: dptdz_threshold !:
209	REAL(wp) :: fac !:
210	REAL(wp) :: height !:
211	REAL(wp) :: pts !:
212	REAL(wp) :: sums_l_eper !:
213	REAL(wp) :: sums_l_etot !:
214	REAL(wp) :: ust !:
215	REAL(wp) :: ust2 !:
216	REAL(wp) :: u2 !:
217	REAL(wp) :: vst !:
218	REAL(wp) :: vst2 !:
219	REAL(wp) :: v2 !:
220	REAL(wp) :: w2 !:
221	REAL(wp) :: z_i(2) !:
222
223	REAL(wp) :: dptdz(nzb+1:nzt+1) !:
224	REAL(wp) :: sums_ll(nzb:nzt+1,2) !:
225
226	CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
227
228	!$acc update host( km, kh, e, pt, qs, qsws, rif, shf, ts, u, usws, v, vsws, w )
229
230	!
231	!-- To be on the safe side, check whether flow_statistics has already been
232	!-- called once after the current time step
233	IF ( flow_statistics_called ) THEN
234
235	message_string = 'flow_statistics is called two times within one ' // &
236	'timestep'
237	CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
238
239	ENDIF
240
241	!
242	!-- Compute statistics for each (sub-)region
243	DO sr = 0, statistic_regions
244
245	!
246	!-- Initialize (local) summation array
247	sums_l = 0.0_wp
248
249	!
250	!-- Store sums that have been computed in other subroutines in summation
251	!-- array
252	sums_l(:,11,:) = sums_l_l(:,sr,:) ! mixing length from diffusivities
253	!-- WARNING: next line still has to be adjusted for OpenMP
254	sums_l(:,21,0) = sums_wsts_bc_l(:,sr) ! heat flux from advec_s_bc
255	sums_l(nzb+9,pr_palm,0) = sums_divold_l(sr) ! old divergence from pres
256	sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr) ! new divergence from pres
257
258	!
259	!-- When calcuating horizontally-averaged total (resolved- plus subgrid-
260	!-- scale) vertical fluxes and velocity variances by using commonly-
261	!-- applied Reynolds-based methods ( e.g. <w'pt'> = (w-<w>)*(pt-<pt>) )
262	!-- in combination with the 5th order advection scheme, pronounced
263	!-- artificial kinks could be observed in the vertical profiles near the
264	!-- surface. Please note: these kinks were not related to the model truth,
265	!-- i.e. these kinks are just related to an evaluation problem.
266	!-- In order avoid these kinks, vertical fluxes and horizontal as well
267	!-- vertical velocity variances are calculated directly within the advection
268	!-- routines, according to the numerical discretization, to evaluate the
269	!-- statistical quantities as they will appear within the prognostic
270	!-- equations.
271	!-- Copy the turbulent quantities, evaluated in the advection routines to
272	!-- the local array sums_l() for further computations.
273	IF ( ws_scheme_mom .AND. sr == 0 ) THEN
274
275	!
276	!-- According to the Neumann bc for the horizontal velocity components,
277	!-- the corresponding fluxes has to satisfiy the same bc.
278	IF ( ocean ) THEN
279	sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
280	sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)
281	ENDIF
282
283	DO i = 0, threads_per_task-1
284	!
285	!-- Swap the turbulent quantities evaluated in advec_ws.
286	sums_l(:,13,i) = sums_wsus_ws_l(:,i) ! wu
287	sums_l(:,15,i) = sums_wsvs_ws_l(:,i) ! wv
288	sums_l(:,30,i) = sums_us2_ws_l(:,i) ! u*2
289	sums_l(:,31,i) = sums_vs2_ws_l(:,i) ! v*2
290	sums_l(:,32,i) = sums_ws2_ws_l(:,i) ! w*2
291	sums_l(:,34,i) = sums_l(:,34,i) + 0.5_wp * &
292	( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) + &
293	sums_ws2_ws_l(:,i) ) ! e*
294	DO k = nzb, nzt
295	sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5_wp * ( &
296	sums_us2_ws_l(k,i) + &
297	sums_vs2_ws_l(k,i) + &
298	sums_ws2_ws_l(k,i) )
299	ENDDO
300	ENDDO
301
302	ENDIF
303
304	IF ( ws_scheme_sca .AND. sr == 0 ) THEN
305
306	DO i = 0, threads_per_task-1
307	sums_l(:,17,i) = sums_wspts_ws_l(:,i) ! wpt from advec_s_ws
308	IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! wsa
309	IF ( humidity .OR. passive_scalar ) sums_l(:,49,i) = &
310	sums_wsqs_ws_l(:,i) !wq
311	ENDDO
312
313	ENDIF
314	!
315	!-- Horizontally averaged profiles of horizontal velocities and temperature.
316	!-- They must have been computed before, because they are already required
317	!-- for other horizontal averages.
318	tn = 0
319
320	!$OMP PARALLEL PRIVATE( i, j, k, tn )
321	#if defined( __intel_openmp_bug )
322	tn = omp_get_thread_num()
323	#else
324	!$ tn = omp_get_thread_num()
325	#endif
326
327	!$OMP DO
328	DO i = nxl, nxr
329	DO j = nys, nyn
330	DO k = nzb_s_inner(j,i), nzt+1
331	sums_l(k,1,tn) = sums_l(k,1,tn) + u(k,j,i) * rmask(j,i,sr)
332	sums_l(k,2,tn) = sums_l(k,2,tn) + v(k,j,i) * rmask(j,i,sr)
333	sums_l(k,4,tn) = sums_l(k,4,tn) + pt(k,j,i) * rmask(j,i,sr)
334	ENDDO
335	ENDDO
336	ENDDO
337
338	!
339	!-- Horizontally averaged profile of salinity
340	IF ( ocean ) THEN
341	!$OMP DO
342	DO i = nxl, nxr
343	DO j = nys, nyn
344	DO k = nzb_s_inner(j,i), nzt+1
345	sums_l(k,23,tn) = sums_l(k,23,tn) + &
346	sa(k,j,i) * rmask(j,i,sr)
347	ENDDO
348	ENDDO
349	ENDDO
350	ENDIF
351
352	!
353	!-- Horizontally averaged profiles of virtual potential temperature,
354	!-- total water content, specific humidity and liquid water potential
355	!-- temperature
356	IF ( humidity ) THEN
357	!$OMP DO
358	DO i = nxl, nxr
359	DO j = nys, nyn
360	DO k = nzb_s_inner(j,i), nzt+1
361	sums_l(k,44,tn) = sums_l(k,44,tn) + &
362	vpt(k,j,i) * rmask(j,i,sr)
363	sums_l(k,41,tn) = sums_l(k,41,tn) + &
364	q(k,j,i) * rmask(j,i,sr)
365	ENDDO
366	ENDDO
367	ENDDO
368	IF ( cloud_physics ) THEN
369	!$OMP DO
370	DO i = nxl, nxr
371	DO j = nys, nyn
372	DO k = nzb_s_inner(j,i), nzt+1
373	sums_l(k,42,tn) = sums_l(k,42,tn) + &
374	( q(k,j,i) - ql(k,j,i) ) * rmask(j,i,sr)
375	sums_l(k,43,tn) = sums_l(k,43,tn) + ( &
376	pt(k,j,i) + l_d_cppt_d_t(k) ql(k,j,i) &
377	) * rmask(j,i,sr)
378	ENDDO
379	ENDDO
380	ENDDO
381	ENDIF
382	ENDIF
383
384	!
385	!-- Horizontally averaged profiles of passive scalar
386	IF ( passive_scalar ) THEN
387	!$OMP DO
388	DO i = nxl, nxr
389	DO j = nys, nyn
390	DO k = nzb_s_inner(j,i), nzt+1
391	sums_l(k,41,tn) = sums_l(k,41,tn) + q(k,j,i) * rmask(j,i,sr)
392	ENDDO
393	ENDDO
394	ENDDO
395	ENDIF
396	!$OMP END PARALLEL
397	!
398	!-- Summation of thread sums
399	IF ( threads_per_task > 1 ) THEN
400	DO i = 1, threads_per_task-1
401	sums_l(:,1,0) = sums_l(:,1,0) + sums_l(:,1,i)
402	sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i)
403	sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i)
404	IF ( ocean ) THEN
405	sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)
406	ENDIF
407	IF ( humidity ) THEN
408	sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
409	sums_l(:,44,0) = sums_l(:,44,0) + sums_l(:,44,i)
410	IF ( cloud_physics ) THEN
411	sums_l(:,42,0) = sums_l(:,42,0) + sums_l(:,42,i)
412	sums_l(:,43,0) = sums_l(:,43,0) + sums_l(:,43,i)
413	ENDIF
414	ENDIF
415	IF ( passive_scalar ) THEN
416	sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
417	ENDIF
418	ENDDO
419	ENDIF
420
421	#if defined( __parallel )
422	!
423	!-- Compute total sum from local sums
424	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
425	CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL, &
426	MPI_SUM, comm2d, ierr )
427	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
428	CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL, &
429	MPI_SUM, comm2d, ierr )
430	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
431	CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL, &
432	MPI_SUM, comm2d, ierr )
433	IF ( ocean ) THEN
434	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
435	CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb, &
436	MPI_REAL, MPI_SUM, comm2d, ierr )
437	ENDIF
438	IF ( humidity ) THEN
439	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
440	CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb, &
441	MPI_REAL, MPI_SUM, comm2d, ierr )
442	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
443	CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
444	MPI_REAL, MPI_SUM, comm2d, ierr )
445	IF ( cloud_physics ) THEN
446	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
447	CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb, &
448	MPI_REAL, MPI_SUM, comm2d, ierr )
449	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
450	CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb, &
451	MPI_REAL, MPI_SUM, comm2d, ierr )
452	ENDIF
453	ENDIF
454
455	IF ( passive_scalar ) THEN
456	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
457	CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
458	MPI_REAL, MPI_SUM, comm2d, ierr )
459	ENDIF
460	#else
461	sums(:,1) = sums_l(:,1,0)
462	sums(:,2) = sums_l(:,2,0)
463	sums(:,4) = sums_l(:,4,0)
464	IF ( ocean ) sums(:,23) = sums_l(:,23,0)
465	IF ( humidity ) THEN
466	sums(:,44) = sums_l(:,44,0)
467	sums(:,41) = sums_l(:,41,0)
468	IF ( cloud_physics ) THEN
469	sums(:,42) = sums_l(:,42,0)
470	sums(:,43) = sums_l(:,43,0)
471	ENDIF
472	ENDIF
473	IF ( passive_scalar ) sums(:,41) = sums_l(:,41,0)
474	#endif
475
476	!
477	!-- Final values are obtained by division by the total number of grid points
478	!-- used for summation. After that store profiles.
479	sums(:,1) = sums(:,1) / ngp_2dh(sr)
480	sums(:,2) = sums(:,2) / ngp_2dh(sr)
481	sums(:,4) = sums(:,4) / ngp_2dh_s_inner(:,sr)
482	hom(:,1,1,sr) = sums(:,1) ! u
483	hom(:,1,2,sr) = sums(:,2) ! v
484	hom(:,1,4,sr) = sums(:,4) ! pt
485
486
487	!
488	!-- Salinity
489	IF ( ocean ) THEN
490	sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr)
491	hom(:,1,23,sr) = sums(:,23) ! sa
492	ENDIF
493
494	!
495	!-- Humidity and cloud parameters
496	IF ( humidity ) THEN
497	sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr)
498	sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
499	hom(:,1,44,sr) = sums(:,44) ! vpt
500	hom(:,1,41,sr) = sums(:,41) ! qv (q)
501	IF ( cloud_physics ) THEN
502	sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr)
503	sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr)
504	hom(:,1,42,sr) = sums(:,42) ! qv
505	hom(:,1,43,sr) = sums(:,43) ! pt
506	ENDIF
507	ENDIF
508
509	!
510	!-- Passive scalar
511	IF ( passive_scalar ) hom(:,1,41,sr) = sums(:,41) / &
512	ngp_2dh_s_inner(:,sr) ! s (q)
513
514	!
515	!-- Horizontally averaged profiles of the remaining prognostic variables,
516	!-- variances, the total and the perturbation energy (single values in last
517	!-- column of sums_l) and some diagnostic quantities.
518	!-- NOTE: for simplicity, nzb_s_inner is used below, although strictly
519	!-- ---- speaking the following k-loop would have to be split up and
520	!-- rearranged according to the staggered grid.
521	!-- However, this implies no error since staggered velocity components
522	!-- are zero at the walls and inside buildings.
523	tn = 0
524	#if defined( __intel_openmp_bug )
525	!$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, &
526	!$OMP tn, ust, ust2, u2, vst, vst2, v2, w2 )
527	tn = omp_get_thread_num()
528	#else
529	!$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2, w2 )
530	!$ tn = omp_get_thread_num()
531	#endif
532	!$OMP DO
533	DO i = nxl, nxr
534	DO j = nys, nyn
535	sums_l_etot = 0.0_wp
536	DO k = nzb_s_inner(j,i), nzt+1
537	!
538	!-- Prognostic and diagnostic variables
539	sums_l(k,3,tn) = sums_l(k,3,tn) + w(k,j,i) * rmask(j,i,sr)
540	sums_l(k,8,tn) = sums_l(k,8,tn) + e(k,j,i) * rmask(j,i,sr)
541	sums_l(k,9,tn) = sums_l(k,9,tn) + km(k,j,i) * rmask(j,i,sr)
542	sums_l(k,10,tn) = sums_l(k,10,tn) + kh(k,j,i) * rmask(j,i,sr)
543	sums_l(k,40,tn) = sums_l(k,40,tn) + p(k,j,i)
544
545	sums_l(k,33,tn) = sums_l(k,33,tn) + &
546	( pt(k,j,i)-hom(k,1,4,sr) )*2 rmask(j,i,sr)
547
548	IF ( humidity ) THEN
549	sums_l(k,70,tn) = sums_l(k,70,tn) + &
550	( q(k,j,i)-hom(k,1,41,sr) )*2 rmask(j,i,sr)
551	ENDIF
552
553	!
554	!-- Higher moments
555	!-- (Computation of the skewness of w further below)
556	sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i)*3 rmask(j,i,sr)
557
558	sums_l_etot = sums_l_etot + &
559	0.5_wp * ( u(k,j,i)2 + v(k,j,i)2 + &
560	w(k,j,i)*2 ) rmask(j,i,sr)
561	ENDDO
562	!
563	!-- Total and perturbation energy for the total domain (being
564	!-- collected in the last column of sums_l). Summation of these
565	!-- quantities is seperated from the previous loop in order to
566	!-- allow vectorization of that loop.
567	sums_l(nzb+4,pr_palm,tn) = sums_l(nzb+4,pr_palm,tn) + sums_l_etot
568	!
569	!-- 2D-arrays (being collected in the last column of sums_l)
570	sums_l(nzb,pr_palm,tn) = sums_l(nzb,pr_palm,tn) + &
571	us(j,i) * rmask(j,i,sr)
572	sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) + &
573	usws(j,i) * rmask(j,i,sr)
574	sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) + &
575	vsws(j,i) * rmask(j,i,sr)
576	sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) + &
577	ts(j,i) * rmask(j,i,sr)
578	IF ( humidity ) THEN
579	sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) + &
580	qs(j,i) * rmask(j,i,sr)
581	ENDIF
582	ENDDO
583	ENDDO
584
585	!
586	!-- Computation of statistics when ws-scheme is not used. Else these
587	!-- quantities are evaluated in the advection routines.
588	IF ( .NOT. ws_scheme_mom .OR. sr /= 0 ) THEN
589	!$OMP DO
590	DO i = nxl, nxr
591	DO j = nys, nyn
592	sums_l_eper = 0.0_wp
593	DO k = nzb_s_inner(j,i), nzt+1
594	u2 = u(k,j,i)**2
595	v2 = v(k,j,i)**2
596	w2 = w(k,j,i)**2
597	ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
598	vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
599
600	sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr)
601	sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr)
602	sums_l(k,32,tn) = sums_l(k,32,tn) + w2 * rmask(j,i,sr)
603	!
604	!-- Perturbation energy
605
606	sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5_wp * &
607	( ust2 + vst2 + w2 ) * rmask(j,i,sr)
608	sums_l_eper = sums_l_eper + &
609	0.5_wp * ( ust2+vst2+w2 ) * rmask(j,i,sr)
610
611	ENDDO
612	sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn) &
613	+ sums_l_eper
614	ENDDO
615	ENDDO
616	ENDIF
617
618	!
619	!-- Horizontally averaged profiles of the vertical fluxes
620
621	!$OMP DO
622	DO i = nxl, nxr
623	DO j = nys, nyn
624	!
625	!-- Subgridscale fluxes (without Prandtl layer from k=nzb,
626	!-- oterwise from k=nzb+1)
627	!-- NOTE: for simplicity, nzb_diff_s_inner is used below, although
628	!-- ---- strictly speaking the following k-loop would have to be
629	!-- split up according to the staggered grid.
630	!-- However, this implies no error since staggered velocity
631	!-- components are zero at the walls and inside buildings.
632
633	DO k = nzb_diff_s_inner(j,i)-1, nzt_diff
634	!
635	!-- Momentum flux w"u"
636	sums_l(k,12,tn) = sums_l(k,12,tn) - 0.25_wp * ( &
637	km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
638	) * ( &
639	( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
640	+ ( w(k,j,i) - w(k,j,i-1) ) * ddx &
641	) * rmask(j,i,sr)
642	!
643	!-- Momentum flux w"v"
644	sums_l(k,14,tn) = sums_l(k,14,tn) - 0.25_wp * ( &
645	km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
646	) * ( &
647	( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
648	+ ( w(k,j,i) - w(k,j-1,i) ) * ddy &
649	) * rmask(j,i,sr)
650	!
651	!-- Heat flux w"pt"
652	sums_l(k,16,tn) = sums_l(k,16,tn) &
653	- 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
654	* ( pt(k+1,j,i) - pt(k,j,i) ) &
655	* ddzu(k+1) * rmask(j,i,sr)
656
657
658	!
659	!-- Salinity flux w"sa"
660	IF ( ocean ) THEN
661	sums_l(k,65,tn) = sums_l(k,65,tn) &
662	- 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
663	* ( sa(k+1,j,i) - sa(k,j,i) ) &
664	* ddzu(k+1) * rmask(j,i,sr)
665	ENDIF
666
667	!
668	!-- Buoyancy flux, water flux (humidity flux) w"q"
669	IF ( humidity ) THEN
670	sums_l(k,45,tn) = sums_l(k,45,tn) &
671	- 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
672	* ( vpt(k+1,j,i) - vpt(k,j,i) ) &
673	* ddzu(k+1) * rmask(j,i,sr)
674	sums_l(k,48,tn) = sums_l(k,48,tn) &
675	- 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
676	* ( q(k+1,j,i) - q(k,j,i) ) &
677	* ddzu(k+1) * rmask(j,i,sr)
678
679	IF ( cloud_physics ) THEN
680	sums_l(k,51,tn) = sums_l(k,51,tn) &
681	- 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
682	* ( ( q(k+1,j,i) - ql(k+1,j,i) )&
683	- ( q(k,j,i) - ql(k,j,i) ) ) &
684	* ddzu(k+1) * rmask(j,i,sr)
685	ENDIF
686	ENDIF
687
688	!
689	!-- Passive scalar flux
690	IF ( passive_scalar ) THEN
691	sums_l(k,48,tn) = sums_l(k,48,tn) &
692	- 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
693	* ( q(k+1,j,i) - q(k,j,i) ) &
694	* ddzu(k+1) * rmask(j,i,sr)
695	ENDIF
696
697	ENDDO
698
699	!
700	!-- Subgridscale fluxes in the Prandtl layer
701	IF ( use_surface_fluxes ) THEN
702	sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
703	usws(j,i) * rmask(j,i,sr) ! w"u"
704	sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
705	vsws(j,i) * rmask(j,i,sr) ! w"v"
706	sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
707	shf(j,i) * rmask(j,i,sr) ! w"pt"
708	sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
709	0.0_wp * rmask(j,i,sr) ! u"pt"
710	sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
711	0.0_wp * rmask(j,i,sr) ! v"pt"
712	IF ( ocean ) THEN
713	sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
714	saswsb(j,i) * rmask(j,i,sr) ! w"sa"
715	ENDIF
716	IF ( humidity ) THEN
717	sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
718	qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
719	sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + ( &
720	( 1.0_wp + 0.61_wp * q(nzb,j,i) ) * &
721	shf(j,i) + 0.61_wp * pt(nzb,j,i) * &
722	qsws(j,i) )
723	IF ( cloud_droplets ) THEN
724	sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + ( &
725	( 1.0_wp + 0.61_wp * q(nzb,j,i) - &
726	ql(nzb,j,i) ) * shf(j,i) + &
727	0.61_wp * pt(nzb,j,i) * qsws(j,i) )
728	ENDIF
729	IF ( cloud_physics ) THEN
730	!
731	!-- Formula does not work if ql(nzb) /= 0.0
732	sums_l(nzb,51,tn) = sums_l(nzb,51,tn) + & ! w"q" (w"qv")
733	qsws(j,i) * rmask(j,i,sr)
734	ENDIF
735	ENDIF
736	IF ( passive_scalar ) THEN
737	sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
738	qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
739	ENDIF
740	ENDIF
741
742	IF ( land_surface ) THEN
743	sums_l(nzb,93,tn) = sums_l(nzb,93,tn) + ghf_eb(j,i)
744	sums_l(nzb,94,tn) = sums_l(nzb,94,tn) + shf_eb(j,i)
745	sums_l(nzb,95,tn) = sums_l(nzb,95,tn) + qsws_eb(j,i)
746	sums_l(nzb,96,tn) = sums_l(nzb,96,tn) + qsws_liq_eb(j,i)
747	sums_l(nzb,97,tn) = sums_l(nzb,97,tn) + qsws_soil_eb(j,i)
748	sums_l(nzb,98,tn) = sums_l(nzb,98,tn) + qsws_veg_eb(j,i)
749	sums_l(nzb,99,tn) = sums_l(nzb,99,tn) + r_a(j,i)
750	sums_l(nzb,100,tn) = sums_l(nzb,100,tn)+ r_s(j,i)
751	ENDIF
752
753	IF ( radiation ) THEN
754	sums_l(nzb,101,tn) = sums_l(nzb,101,tn) + rad_net(j,i)
755	sums_l(nzb,102,tn) = sums_l(nzb,102,tn) + rad_lw_in(nzb,j,i)
756	sums_l(nzb,103,tn) = sums_l(nzb,103,tn) + rad_lw_out(nzb,j,i)
757	sums_l(nzb,104,tn) = sums_l(nzb,104,tn) + rad_sw_in(nzb,j,i)
758	sums_l(nzb,105,tn) = sums_l(nzb,105,tn) + rad_sw_out(nzb,j,i)
759
760	#if defined ( __rrtmg )
761	IF ( radiation_scheme == 'rrtmg' ) THEN
762	sums_l(nzb,106,tn) = sums_l(nzb,106,tn) + rrtm_aldif(0,j,i)
763	sums_l(nzb,107,tn) = sums_l(nzb,107,tn) + rrtm_aldir(0,j,i)
764	sums_l(nzb,108,tn) = sums_l(nzb,108,tn) + rrtm_asdif(0,j,i)
765	sums_l(nzb,109,tn) = sums_l(nzb,109,tn) + rrtm_asdir(0,j,i)
766	ENDIF
767	#endif
768
769	ENDIF
770
771	!
772	!-- Subgridscale fluxes at the top surface
773	IF ( use_top_fluxes ) THEN
774	sums_l(nzt:nzt+1,12,tn) = sums_l(nzt:nzt+1,12,tn) + &
775	uswst(j,i) * rmask(j,i,sr) ! w"u"
776	sums_l(nzt:nzt+1,14,tn) = sums_l(nzt:nzt+1,14,tn) + &
777	vswst(j,i) * rmask(j,i,sr) ! w"v"
778	sums_l(nzt:nzt+1,16,tn) = sums_l(nzt:nzt+1,16,tn) + &
779	tswst(j,i) * rmask(j,i,sr) ! w"pt"
780	sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + &
781	0.0_wp * rmask(j,i,sr) ! u"pt"
782	sums_l(nzt:nzt+1,61,tn) = sums_l(nzt:nzt+1,61,tn) + &
783	0.0_wp * rmask(j,i,sr) ! v"pt"
784
785	IF ( ocean ) THEN
786	sums_l(nzt,65,tn) = sums_l(nzt,65,tn) + &
787	saswst(j,i) * rmask(j,i,sr) ! w"sa"
788	ENDIF
789	IF ( humidity ) THEN
790	sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
791	qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
792	sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + ( &
793	( 1.0_wp + 0.61_wp * q(nzt,j,i) ) * &
794	tswst(j,i) + 0.61_wp * pt(nzt,j,i) * &
795	qswst(j,i) )
796	IF ( cloud_droplets ) THEN
797	sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + ( &
798	( 1.0_wp + 0.61_wp * q(nzt,j,i) - &
799	ql(nzt,j,i) ) * tswst(j,i) + &
800	0.61_wp * pt(nzt,j,i) * qswst(j,i) )
801	ENDIF
802	IF ( cloud_physics ) THEN
803	!
804	!-- Formula does not work if ql(nzb) /= 0.0
805	sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + & ! w"q" (w"qv")
806	qswst(j,i) * rmask(j,i,sr)
807	ENDIF
808	ENDIF
809	IF ( passive_scalar ) THEN
810	sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
811	qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
812	ENDIF
813	ENDIF
814
815	!
816	!-- Resolved fluxes (can be computed for all horizontal points)
817	!-- NOTE: for simplicity, nzb_s_inner is used below, although strictly
818	!-- ---- speaking the following k-loop would have to be split up and
819	!-- rearranged according to the staggered grid.
820	DO k = nzb_s_inner(j,i), nzt
821	ust = 0.5_wp * ( u(k,j,i) - hom(k,1,1,sr) + &
822	u(k+1,j,i) - hom(k+1,1,1,sr) )
823	vst = 0.5_wp * ( v(k,j,i) - hom(k,1,2,sr) + &
824	v(k+1,j,i) - hom(k+1,1,2,sr) )
825	pts = 0.5_wp * ( pt(k,j,i) - hom(k,1,4,sr) + &
826	pt(k+1,j,i) - hom(k+1,1,4,sr) )
827
828	!-- Higher moments
829	sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)*2 &
830	rmask(j,i,sr)
831	sums_l(k,36,tn) = sums_l(k,36,tn) + pts*2 w(k,j,i) * &
832	rmask(j,i,sr)
833
834	!
835	!-- Salinity flux and density (density does not belong to here,
836	!-- but so far there is no other suitable place to calculate)
837	IF ( ocean ) THEN
838	IF( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN
839	pts = 0.5_wp * ( sa(k,j,i) - hom(k,1,23,sr) + &
840	sa(k+1,j,i) - hom(k+1,1,23,sr) )
841	sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) * &
842	rmask(j,i,sr)
843	ENDIF
844	sums_l(k,64,tn) = sums_l(k,64,tn) + rho(k,j,i) * &
845	rmask(j,i,sr)
846	sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) * &
847	rmask(j,i,sr)
848	ENDIF
849
850	!
851	!-- Buoyancy flux, water flux, humidity flux, liquid water
852	!-- content, rain drop concentration and rain water content
853	IF ( humidity ) THEN
854	IF ( cloud_physics .OR. cloud_droplets ) THEN
855	pts = 0.5_wp * ( vpt(k,j,i) - hom(k,1,44,sr) + &
856	vpt(k+1,j,i) - hom(k+1,1,44,sr) )
857	sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
858	rmask(j,i,sr)
859	IF ( .NOT. cloud_droplets ) THEN
860	pts = 0.5_wp * &
861	( ( q(k,j,i) - ql(k,j,i) ) - &
862	hom(k,1,42,sr) + &
863	( q(k+1,j,i) - ql(k+1,j,i) ) - &
864	hom(k+1,1,42,sr) )
865	sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) * &
866	rmask(j,i,sr)
867	IF ( icloud_scheme == 0 ) THEN
868	sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * &
869	rmask(j,i,sr)
870	sums_l(k,75,tn) = sums_l(k,75,tn) + qc(k,j,i) * &
871	rmask(j,i,sr)
872	IF ( precipitation ) THEN
873	sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) * &
874	rmask(j,i,sr)
875	sums_l(k,74,tn) = sums_l(k,74,tn) + qr(k,j,i) * &
876	rmask(j,i,sr)
877	sums_l(k,76,tn) = sums_l(k,76,tn) + prr(k,j,i) *&
878	rmask(j,i,sr)
879	ENDIF
880	ELSE
881	sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * &
882	rmask(j,i,sr)
883	ENDIF
884	ELSE
885	sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * &
886	rmask(j,i,sr)
887	ENDIF
888	ELSE
889	IF( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN
890	pts = 0.5_wp * ( vpt(k,j,i) - hom(k,1,44,sr) + &
891	vpt(k+1,j,i) - hom(k+1,1,44,sr) )
892	sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
893	rmask(j,i,sr)
894	ELSE IF ( ws_scheme_sca .AND. sr == 0 ) THEN
895	sums_l(k,46,tn) = ( 1.0_wp + 0.61_wp * &
896	hom(k,1,41,sr) ) * &
897	sums_l(k,17,tn) + &
898	0.61_wp * hom(k,1,4,sr) * &
899	sums_l(k,49,tn)
900	END IF
901	END IF
902	ENDIF
903	!
904	!-- Passive scalar flux
905	IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca &
906	.OR. sr /= 0 ) ) THEN
907	pts = 0.5_wp * ( q(k,j,i) - hom(k,1,41,sr) + &
908	q(k+1,j,i) - hom(k+1,1,41,sr) )
909	sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
910	rmask(j,i,sr)
911	ENDIF
912
913	!
914	!-- Energy flux we
915	!-- has to be adjusted
916	sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5_wp * &
917	( ust2 + vst2 + w(k,j,i)**2 ) &
918	* rmask(j,i,sr)
919	ENDDO
920	ENDDO
921	ENDDO
922	!
923	!-- For speed optimization fluxes which have been computed in part directly
924	!-- inside the WS advection routines are treated seperatly
925	!-- Momentum fluxes first:
926	IF ( .NOT. ws_scheme_mom .OR. sr /= 0 ) THEN
927	!$OMP DO
928	DO i = nxl, nxr
929	DO j = nys, nyn
930	DO k = nzb_diff_s_inner(j,i)-1, nzt_diff
931	ust = 0.5_wp * ( u(k,j,i) - hom(k,1,1,sr) + &
932	u(k+1,j,i) - hom(k+1,1,1,sr) )
933	vst = 0.5_wp * ( v(k,j,i) - hom(k,1,2,sr) + &
934	v(k+1,j,i) - hom(k+1,1,2,sr) )
935	!
936	!-- Momentum flux wu
937	sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5_wp * &
938	( w(k,j,i-1) + w(k,j,i) ) &
939	* ust * rmask(j,i,sr)
940	!
941	!-- Momentum flux wv
942	sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5_wp * &
943	( w(k,j-1,i) + w(k,j,i) ) &
944	* vst * rmask(j,i,sr)
945	ENDDO
946	ENDDO
947	ENDDO
948
949	ENDIF
950	IF ( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN
951	!$OMP DO
952	DO i = nxl, nxr
953	DO j = nys, nyn
954	DO k = nzb_diff_s_inner(j,i)-1, nzt_diff
955	!
956	!-- Vertical heat flux
957	sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5_wp * &
958	( pt(k,j,i) - hom(k,1,4,sr) + &
959	pt(k+1,j,i) - hom(k+1,1,4,sr) ) &
960	* w(k,j,i) * rmask(j,i,sr)
961	IF ( humidity ) THEN
962	pts = 0.5_wp * ( q(k,j,i) - hom(k,1,41,sr) + &
963	q(k+1,j,i) - hom(k+1,1,41,sr) )
964	sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
965	rmask(j,i,sr)
966	ENDIF
967	ENDDO
968	ENDDO
969	ENDDO
970
971	ENDIF
972
973	!
974	!-- Density at top follows Neumann condition
975	IF ( ocean ) THEN
976	sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
977	sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
978	ENDIF
979
980	!
981	!-- Divergence of vertical flux of resolved scale energy and pressure
982	!-- fluctuations as well as flux of pressure fluctuation itself (68).
983	!-- First calculate the products, then the divergence.
984	!-- Calculation is time consuming. Do it only, if profiles shall be plotted.
985	IF ( hom(nzb+1,2,55,0) /= 0.0_wp .OR. hom(nzb+1,2,68,0) /= 0.0_wp ) THEN
986
987	sums_ll = 0.0_wp ! local array
988
989	!$OMP DO
990	DO i = nxl, nxr
991	DO j = nys, nyn
992	DO k = nzb_s_inner(j,i)+1, nzt
993
994	sums_ll(k,1) = sums_ll(k,1) + 0.5_wp * w(k,j,i) * ( &
995	( 0.25_wp * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1) ) &
996	- 0.5_wp * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) ) )**2&
997	+ ( 0.25_wp * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i) ) &
998	- 0.5_wp * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) ) )**2& &
999	+ w(k,j,i)**2 )
1000
1001	sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i) &
1002	* ( p(k,j,i) + p(k+1,j,i) )
1003
1004	ENDDO
1005	ENDDO
1006	ENDDO
1007	sums_ll(0,1) = 0.0_wp ! because w is zero at the bottom
1008	sums_ll(nzt+1,1) = 0.0_wp
1009	sums_ll(0,2) = 0.0_wp
1010	sums_ll(nzt+1,2) = 0.0_wp
1011
1012	DO k = nzb+1, nzt
1013	sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
1014	sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
1015	sums_l(k,68,tn) = sums_ll(k,2)
1016	ENDDO
1017	sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
1018	sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
1019	sums_l(nzb,68,tn) = 0.0_wp ! because w* = 0 at nzb
1020
1021	ENDIF
1022
1023	!
1024	!-- Divergence of vertical flux of SGS TKE and the flux itself (69)
1025	IF ( hom(nzb+1,2,57,0) /= 0.0_wp .OR. hom(nzb+1,2,69,0) /= 0.0_wp ) THEN
1026
1027	!$OMP DO
1028	DO i = nxl, nxr
1029	DO j = nys, nyn
1030	DO k = nzb_s_inner(j,i)+1, nzt
1031
1032	sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5_wp * ( &
1033	(km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
1034	- (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k) &
1035	) * ddzw(k)
1036
1037	sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5_wp * ( &
1038	(km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
1039	)
1040
1041	ENDDO
1042	ENDDO
1043	ENDDO
1044	sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
1045	sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
1046
1047	ENDIF
1048
1049	!
1050	!-- Horizontal heat fluxes (subgrid, resolved, total).
1051	!-- Do it only, if profiles shall be plotted.
1052	IF ( hom(nzb+1,2,58,0) /= 0.0_wp ) THEN
1053
1054	!$OMP DO
1055	DO i = nxl, nxr
1056	DO j = nys, nyn
1057	DO k = nzb_s_inner(j,i)+1, nzt
1058	!
1059	!-- Subgrid horizontal heat fluxes u"pt", v"pt"
1060	sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5_wp * &
1061	( kh(k,j,i) + kh(k,j,i-1) ) &
1062	* ( pt(k,j,i-1) - pt(k,j,i) ) &
1063	* ddx * rmask(j,i,sr)
1064	sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp * &
1065	( kh(k,j,i) + kh(k,j-1,i) ) &
1066	* ( pt(k,j-1,i) - pt(k,j,i) ) &
1067	* ddy * rmask(j,i,sr)
1068	!
1069	!-- Resolved horizontal heat fluxes upt, vpt
1070	sums_l(k,59,tn) = sums_l(k,59,tn) + &
1071	( u(k,j,i) - hom(k,1,1,sr) ) &
1072	* 0.5_wp * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
1073	pt(k,j,i) - hom(k,1,4,sr) )
1074	pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
1075	pt(k,j,i) - hom(k,1,4,sr) )
1076	sums_l(k,62,tn) = sums_l(k,62,tn) + &
1077	( v(k,j,i) - hom(k,1,2,sr) ) &
1078	* 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
1079	pt(k,j,i) - hom(k,1,4,sr) )
1080	ENDDO
1081	ENDDO
1082	ENDDO
1083	!
1084	!-- Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
1085	sums_l(nzb,58,tn) = 0.0_wp
1086	sums_l(nzb,59,tn) = 0.0_wp
1087	sums_l(nzb,60,tn) = 0.0_wp
1088	sums_l(nzb,61,tn) = 0.0_wp
1089	sums_l(nzb,62,tn) = 0.0_wp
1090	sums_l(nzb,63,tn) = 0.0_wp
1091
1092	ENDIF
1093
1094	!
1095	!-- Collect current large scale advection and subsidence tendencies for
1096	!-- data output
1097	IF ( large_scale_forcing .AND. ( simulated_time .GT. 0.0_wp ) ) THEN
1098	!
1099	!-- Interpolation in time of LSF_DATA
1100	nt = 1
1101	DO WHILE ( simulated_time - dt_3d > time_vert(nt) )
1102	nt = nt + 1
1103	ENDDO
1104	IF ( simulated_time - dt_3d /= time_vert(nt) ) THEN
1105	nt = nt - 1
1106	ENDIF
1107
1108	fac = ( simulated_time - dt_3d - time_vert(nt) ) &
1109	/ ( time_vert(nt+1)-time_vert(nt) )
1110
1111
1112	DO k = nzb, nzt
1113	sums_ls_l(k,0) = td_lsa_lpt(k,nt) &
1114	+ fac * ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )
1115	sums_ls_l(k,1) = td_lsa_q(k,nt) &
1116	+ fac * ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) )
1117	ENDDO
1118
1119	sums_ls_l(nzt+1,0) = sums_ls_l(nzt,0)
1120	sums_ls_l(nzt+1,1) = sums_ls_l(nzt,1)
1121
1122	IF ( large_scale_subsidence .AND. use_subsidence_tendencies ) THEN
1123
1124	DO k = nzb, nzt
1125	sums_ls_l(k,2) = td_sub_lpt(k,nt) + fac * &
1126	( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )
1127	sums_ls_l(k,3) = td_sub_q(k,nt) + fac * &
1128	( td_sub_q(k,nt+1) - td_sub_q(k,nt) )
1129	ENDDO
1130
1131	sums_ls_l(nzt+1,2) = sums_ls_l(nzt,2)
1132	sums_ls_l(nzt+1,3) = sums_ls_l(nzt,3)
1133
1134	ENDIF
1135
1136	ENDIF
1137
1138
1139	IF ( land_surface ) THEN
1140	!$OMP DO
1141	DO i = nxl, nxr
1142	DO j = nys, nyn
1143	DO k = nzb_soil, nzt_soil
1144	sums_l(k,89,tn) = sums_l(k,89,tn) + t_soil(k,j,i) * rmask(j,i,sr)
1145	sums_l(k,91,tn) = sums_l(k,91,tn) + m_soil(k,j,i) * rmask(j,i,sr)
1146	ENDDO
1147	ENDDO
1148	ENDDO
1149	ENDIF
1150
1151	IF ( radiation .AND. radiation_scheme == 'rrtmg' ) THEN
1152	!$OMP DO
1153	DO i = nxl, nxr
1154	DO j = nys, nyn
1155	DO k = nzb_s_inner(j,i)+1, nzt+1
1156	sums_l(k,102,tn) = sums_l(k,102,tn) + rad_lw_in(k,j,i) * rmask(j,i,sr)
1157	sums_l(k,103,tn) = sums_l(k,103,tn) + rad_lw_out(k,j,i) * rmask(j,i,sr)
1158	sums_l(k,104,tn) = sums_l(k,104,tn) + rad_sw_in(k,j,i) * rmask(j,i,sr)
1159	sums_l(k,105,tn) = sums_l(k,105,tn) + rad_sw_out(k,j,i) * rmask(j,i,sr)
1160	ENDDO
1161	ENDDO
1162	ENDDO
1163	ENDIF
1164	!
1165	!-- Calculate the user-defined profiles
1166	CALL user_statistics( 'profiles', sr, tn )
1167	!$OMP END PARALLEL
1168
1169	!
1170	!-- Summation of thread sums
1171	IF ( threads_per_task > 1 ) THEN
1172	DO i = 1, threads_per_task-1
1173	sums_l(:,3,0) = sums_l(:,3,0) + sums_l(:,3,i)
1174	sums_l(:,4:40,0) = sums_l(:,4:40,0) + sums_l(:,4:40,i)
1175	sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
1176	sums_l(:,45:pr_palm,i)
1177	IF ( max_pr_user > 0 ) THEN
1178	sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
1179	sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
1180	sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
1181	ENDIF
1182	ENDDO
1183	ENDIF
1184
1185	#if defined( __parallel )
1186
1187	!
1188	!-- Compute total sum from local sums
1189	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1190	CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, &
1191	MPI_SUM, comm2d, ierr )
1192	IF ( large_scale_forcing ) THEN
1193	CALL MPI_ALLREDUCE( sums_ls_l(nzb,2), sums(nzb,83), ngp_sums_ls, &
1194	MPI_REAL, MPI_SUM, comm2d, ierr )
1195	ENDIF
1196	#else
1197	sums = sums_l(:,:,0)
1198	IF ( large_scale_forcing ) THEN
1199	sums(:,81:88) = sums_ls_l
1200	ENDIF
1201	#endif
1202
1203	!
1204	!-- Final values are obtained by division by the total number of grid points
1205	!-- used for summation. After that store profiles.
1206	!-- Profiles:
1207	DO k = nzb, nzt+1
1208	sums(k,3) = sums(k,3) / ngp_2dh(sr)
1209	sums(k,8:11) = sums(k,8:11) / ngp_2dh_s_inner(k,sr)
1210	sums(k,12:22) = sums(k,12:22) / ngp_2dh(sr)
1211	sums(k,23:29) = sums(k,23:29) / ngp_2dh_s_inner(k,sr)
1212	sums(k,30:32) = sums(k,30:32) / ngp_2dh(sr)
1213	sums(k,33:34) = sums(k,33:34) / ngp_2dh_s_inner(k,sr)
1214	sums(k,35:39) = sums(k,35:39) / ngp_2dh(sr)
1215	sums(k,40) = sums(k,40) / ngp_2dh_s_inner(k,sr)
1216	sums(k,45:53) = sums(k,45:53) / ngp_2dh(sr)
1217	sums(k,54) = sums(k,54) / ngp_2dh_s_inner(k,sr)
1218	sums(k,55:63) = sums(k,55:63) / ngp_2dh(sr)
1219	sums(k,64) = sums(k,64) / ngp_2dh_s_inner(k,sr)
1220	sums(k,65:69) = sums(k,65:69) / ngp_2dh(sr)
1221	sums(k,70:80) = sums(k,70:80) / ngp_2dh_s_inner(k,sr)
1222	sums(k,81:88) = sums(k,81:88) / ngp_2dh(sr)
1223	sums(k,89:105) = sums(k,89:105) / ngp_2dh(sr)
1224	sums(k,106:pr_palm-2) = sums(k,106:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
1225	ENDDO
1226
1227	!-- Upstream-parts
1228	sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
1229	!-- u* and so on
1230	!-- As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
1231	!-- size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
1232	!-- above the topography, they are being divided by ngp_2dh(sr)
1233	sums(nzb:nzb+3,pr_palm) = sums(nzb:nzb+3,pr_palm) / &
1234	ngp_2dh(sr)
1235	sums(nzb+12,pr_palm) = sums(nzb+12,pr_palm) / & ! qs
1236	ngp_2dh(sr)
1237	!-- eges, e*
1238	sums(nzb+4:nzb+5,pr_palm) = sums(nzb+4:nzb+5,pr_palm) / &
1239	ngp_3d(sr)
1240	!-- Old and new divergence
1241	sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
1242	ngp_3d_inner(sr)
1243
1244	!-- User-defined profiles
1245	IF ( max_pr_user > 0 ) THEN
1246	DO k = nzb, nzt+1
1247	sums(k,pr_palm+1:pr_palm+max_pr_user) = &
1248	sums(k,pr_palm+1:pr_palm+max_pr_user) / &
1249	ngp_2dh_s_inner(k,sr)
1250	ENDDO
1251	ENDIF
1252
1253	!
1254	!-- Collect horizontal average in hom.
1255	!-- Compute deduced averages (e.g. total heat flux)
1256	hom(:,1,3,sr) = sums(:,3) ! w
1257	hom(:,1,8,sr) = sums(:,8) ! e profiles 5-7 are initial profiles
1258	hom(:,1,9,sr) = sums(:,9) ! km
1259	hom(:,1,10,sr) = sums(:,10) ! kh
1260	hom(:,1,11,sr) = sums(:,11) ! l
1261	hom(:,1,12,sr) = sums(:,12) ! w"u"
1262	hom(:,1,13,sr) = sums(:,13) ! wu
1263	hom(:,1,14,sr) = sums(:,14) ! w"v"
1264	hom(:,1,15,sr) = sums(:,15) ! wv
1265	hom(:,1,16,sr) = sums(:,16) ! w"pt"
1266	hom(:,1,17,sr) = sums(:,17) ! wpt
1267	hom(:,1,18,sr) = sums(:,16) + sums(:,17) ! wpt
1268	hom(:,1,19,sr) = sums(:,12) + sums(:,13) ! wu
1269	hom(:,1,20,sr) = sums(:,14) + sums(:,15) ! wv
1270	hom(:,1,21,sr) = sums(:,21) ! wptBC
1271	hom(:,1,22,sr) = sums(:,16) + sums(:,21) ! wptBC
1272	! profile 24 is initial profile (sa)
1273	! profiles 25-29 left empty for initial
1274	! profiles
1275	hom(:,1,30,sr) = sums(:,30) ! u*2
1276	hom(:,1,31,sr) = sums(:,31) ! v*2
1277	hom(:,1,32,sr) = sums(:,32) ! w*2
1278	hom(:,1,33,sr) = sums(:,33) ! pt*2
1279	hom(:,1,34,sr) = sums(:,34) ! e*
1280	hom(:,1,35,sr) = sums(:,35) ! w2pt
1281	hom(:,1,36,sr) = sums(:,36) ! wpt2
1282	hom(:,1,37,sr) = sums(:,37) ! we
1283	hom(:,1,38,sr) = sums(:,38) ! w*3
1284	hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20_wp )**1.5_wp ! Sw
1285	hom(:,1,40,sr) = sums(:,40) ! p
1286	hom(:,1,45,sr) = sums(:,45) ! w"vpt"
1287	hom(:,1,46,sr) = sums(:,46) ! wvpt
1288	hom(:,1,47,sr) = sums(:,45) + sums(:,46) ! wvpt
1289	hom(:,1,48,sr) = sums(:,48) ! w"q" (w"qv")
1290	hom(:,1,49,sr) = sums(:,49) ! wq (wqv)
1291	hom(:,1,50,sr) = sums(:,48) + sums(:,49) ! wq (wqv)
1292	hom(:,1,51,sr) = sums(:,51) ! w"qv"
1293	hom(:,1,52,sr) = sums(:,52) ! wqv
1294	hom(:,1,53,sr) = sums(:,52) + sums(:,51) ! wq (wqv)
1295	hom(:,1,54,sr) = sums(:,54) ! ql
1296	hom(:,1,55,sr) = sums(:,55) ! wuu*/dz
1297	hom(:,1,56,sr) = sums(:,56) ! wp/dz
1298	hom(:,1,57,sr) = sums(:,57) ! ( w"e + w"p"/rho )/dz
1299	hom(:,1,58,sr) = sums(:,58) ! u"pt"
1300	hom(:,1,59,sr) = sums(:,59) ! upt
1301	hom(:,1,60,sr) = sums(:,58) + sums(:,59) ! upt_t
1302	hom(:,1,61,sr) = sums(:,61) ! v"pt"
1303	hom(:,1,62,sr) = sums(:,62) ! vpt
1304	hom(:,1,63,sr) = sums(:,61) + sums(:,62) ! vpt_t
1305	hom(:,1,64,sr) = sums(:,64) ! rho
1306	hom(:,1,65,sr) = sums(:,65) ! w"sa"
1307	hom(:,1,66,sr) = sums(:,66) ! wsa
1308	hom(:,1,67,sr) = sums(:,65) + sums(:,66) ! wsa
1309	hom(:,1,68,sr) = sums(:,68) ! wp
1310	hom(:,1,69,sr) = sums(:,69) ! w"e + w"p"/rho
1311	hom(:,1,70,sr) = sums(:,70) ! q*2
1312	hom(:,1,71,sr) = sums(:,71) ! prho
1313	hom(:,1,72,sr) = hyp * 1E-4_wp ! hyp in dbar
1314	hom(:,1,73,sr) = sums(:,73) ! nr
1315	hom(:,1,74,sr) = sums(:,74) ! qr
1316	hom(:,1,75,sr) = sums(:,75) ! qc
1317	hom(:,1,76,sr) = sums(:,76) ! prr (precipitation rate)
1318	! 77 is initial density profile
1319	hom(:,1,78,sr) = ug ! ug
1320	hom(:,1,79,sr) = vg ! vg
1321	hom(:,1,80,sr) = w_subs ! w_subs
1322
1323	IF ( large_scale_forcing ) THEN
1324	hom(:,1,81,sr) = sums_ls_l(:,0) ! td_lsa_lpt
1325	hom(:,1,82,sr) = sums_ls_l(:,1) ! td_lsa_q
1326	IF ( use_subsidence_tendencies ) THEN
1327	hom(:,1,83,sr) = sums_ls_l(:,2) ! td_sub_lpt
1328	hom(:,1,84,sr) = sums_ls_l(:,3) ! td_sub_q
1329	ELSE
1330	hom(:,1,83,sr) = sums(:,83) ! td_sub_lpt
1331	hom(:,1,84,sr) = sums(:,84) ! td_sub_q
1332	ENDIF
1333	hom(:,1,85,sr) = sums(:,85) ! td_nud_lpt
1334	hom(:,1,86,sr) = sums(:,86) ! td_nud_q
1335	hom(:,1,87,sr) = sums(:,87) ! td_nud_u
1336	hom(:,1,88,sr) = sums(:,88) ! td_nud_v
1337	ENDIF
1338
1339	IF ( land_surface ) THEN
1340	hom(:,1,89,sr) = sums(:,89) ! t_soil
1341	! 90 is initial t_soil profile
1342	hom(:,1,91,sr) = sums(:,91) ! m_soil
1343	! 92 is initial m_soil profile
1344	hom(:,1,93,sr) = sums(:,93) ! ghf_eb
1345	hom(:,1,94,sr) = sums(:,94) ! shf_eb
1346	hom(:,1,95,sr) = sums(:,95) ! qsws_eb
1347	hom(:,1,96,sr) = sums(:,96) ! qsws_liq_eb
1348	hom(:,1,97,sr) = sums(:,97) ! qsws_soil_eb
1349	hom(:,1,98,sr) = sums(:,98) ! qsws_veg_eb
1350	hom(:,1,99,sr) = sums(:,99) ! r_a
1351	hom(:,1,100,sr) = sums(:,100) ! r_s
1352
1353	ENDIF
1354
1355	IF ( radiation ) THEN
1356	hom(:,1,101,sr) = sums(:,101) ! rad_net
1357	hom(:,1,102,sr) = sums(:,102) ! rad_lw_in
1358	hom(:,1,103,sr) = sums(:,103) ! rad_lw_out
1359	hom(:,1,104,sr) = sums(:,104) ! rad_sw_in
1360	hom(:,1,105,sr) = sums(:,105) ! rad_sw_out
1361
1362	#if defined ( __rrtmg )
1363	IF ( radiation_scheme == 'rrtmg' ) THEN
1364	hom(:,1,106,sr) = sums(:,106) ! rrtm_aldif
1365	hom(:,1,107,sr) = sums(:,107) ! rrtm_aldir
1366	hom(:,1,108,sr) = sums(:,108) ! rrtm_asdif
1367	hom(:,1,109,sr) = sums(:,109) ! rrtm_asdir
1368	ENDIF
1369	#endif
1370
1371	ENDIF
1372
1373	hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
1374	! upstream-parts u_x, u_y, u_z, v_x,
1375	! v_y, usw. (in last but one profile)
1376	hom(:,1,pr_palm,sr) = sums(:,pr_palm)
1377	! u, w'u', w'v', t (in last profile)
1378
1379	IF ( max_pr_user > 0 ) THEN ! user-defined profiles
1380	hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
1381	sums(:,pr_palm+1:pr_palm+max_pr_user)
1382	ENDIF
1383
1384	!
1385	!-- Determine the boundary layer height using two different schemes.
1386	!-- First scheme: Starting from the Earth's (Ocean's) surface, look for the
1387	!-- first relative minimum (maximum) of the total heat flux.
1388	!-- The corresponding height is assumed as the boundary layer height, if it
1389	!-- is less than 1.5 times the height where the heat flux becomes negative
1390	!-- (positive) for the first time.
1391	z_i(1) = 0.0_wp
1392	first = .TRUE.
1393
1394	IF ( ocean ) THEN
1395	DO k = nzt, nzb+1, -1
1396	IF ( first .AND. hom(k,1,18,sr) < 0.0_wp &
1397	.AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp) THEN
1398	first = .FALSE.
1399	height = zw(k)
1400	ENDIF
1401	IF ( hom(k,1,18,sr) < 0.0_wp .AND. &
1402	abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND. &
1403	hom(k-1,1,18,sr) > hom(k,1,18,sr) ) THEN
1404	IF ( zw(k) < 1.5_wp * height ) THEN
1405	z_i(1) = zw(k)
1406	ELSE
1407	z_i(1) = height
1408	ENDIF
1409	EXIT
1410	ENDIF
1411	ENDDO
1412	ELSE
1413	DO k = nzb, nzt-1
1414	IF ( first .AND. hom(k,1,18,sr) < 0.0_wp &
1415	.AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp ) THEN
1416	first = .FALSE.
1417	height = zw(k)
1418	ENDIF
1419	IF ( hom(k,1,18,sr) < 0.0_wp .AND. &
1420	abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND. &
1421	hom(k+1,1,18,sr) > hom(k,1,18,sr) ) THEN
1422	IF ( zw(k) < 1.5_wp * height ) THEN
1423	z_i(1) = zw(k)
1424	ELSE
1425	z_i(1) = height
1426	ENDIF
1427	EXIT
1428	ENDIF
1429	ENDDO
1430	ENDIF
1431
1432	!
1433	!-- Second scheme: Gradient scheme from Sullivan et al. (1998), modified
1434	!-- by Uhlenbrock(2006). The boundary layer height is the height with the
1435	!-- maximal local temperature gradient: starting from the second (the last
1436	!-- but one) vertical gridpoint, the local gradient must be at least
1437	!-- 0.2K/100m and greater than the next four gradients.
1438	!-- WARNING: The threshold value of 0.2K/100m must be adjusted for the
1439	!-- ocean case!
1440	z_i(2) = 0.0_wp
1441	DO k = nzb+1, nzt+1
1442	dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
1443	ENDDO
1444	dptdz_threshold = 0.2_wp / 100.0_wp
1445
1446	IF ( ocean ) THEN
1447	DO k = nzt+1, nzb+5, -1
1448	IF ( dptdz(k) > dptdz_threshold .AND. &
1449	dptdz(k) > dptdz(k-1) .AND. dptdz(k) > dptdz(k-2) .AND. &
1450	dptdz(k) > dptdz(k-3) .AND. dptdz(k) > dptdz(k-4) ) THEN
1451	z_i(2) = zw(k-1)
1452	EXIT
1453	ENDIF
1454	ENDDO
1455	ELSE
1456	DO k = nzb+1, nzt-3
1457	IF ( dptdz(k) > dptdz_threshold .AND. &
1458	dptdz(k) > dptdz(k+1) .AND. dptdz(k) > dptdz(k+2) .AND. &
1459	dptdz(k) > dptdz(k+3) .AND. dptdz(k) > dptdz(k+4) ) THEN
1460	z_i(2) = zw(k-1)
1461	EXIT
1462	ENDIF
1463	ENDDO
1464	ENDIF
1465
1466	hom(nzb+6,1,pr_palm,sr) = z_i(1)
1467	hom(nzb+7,1,pr_palm,sr) = z_i(2)
1468
1469	!
1470	!-- Computation of both the characteristic vertical velocity and
1471	!-- the characteristic convective boundary layer temperature.
1472	!-- The horizontal average at nzb+1 is input for the average temperature.
1473	IF ( hom(nzb,1,18,sr) > 0.0_wp .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8_wp &
1474	.AND. z_i(1) /= 0.0_wp ) THEN
1475	hom(nzb+8,1,pr_palm,sr) = ( g / hom(nzb+1,1,4,sr) * &
1476	hom(nzb,1,18,sr) * &
1477	ABS( z_i(1) ) )**0.333333333_wp
1478	!-- so far this only works if Prandtl layer is used
1479	hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
1480	ELSE
1481	hom(nzb+8,1,pr_palm,sr) = 0.0_wp
1482	hom(nzb+11,1,pr_palm,sr) = 0.0_wp
1483	ENDIF
1484
1485	!
1486	!-- Collect the time series quantities
1487	ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr) ! E
1488	ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr) ! E*
1489	ts_value(3,sr) = dt_3d
1490	ts_value(4,sr) = hom(nzb,1,pr_palm,sr) ! u*
1491	ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr) ! th*
1492	ts_value(6,sr) = u_max
1493	ts_value(7,sr) = v_max
1494	ts_value(8,sr) = w_max
1495	ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr) ! new divergence
1496	ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr) ! old Divergence
1497	ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr) ! z_i(1)
1498	ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr) ! z_i(2)
1499	ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr) ! w*
1500	ts_value(14,sr) = hom(nzb,1,16,sr) ! w'pt' at k=0
1501	ts_value(15,sr) = hom(nzb+1,1,16,sr) ! w'pt' at k=1
1502	ts_value(16,sr) = hom(nzb+1,1,18,sr) ! wpt at k=1
1503	ts_value(17,sr) = hom(nzb,1,4,sr) ! pt(0)
1504	ts_value(18,sr) = hom(nzb+1,1,4,sr) ! pt(zp)
1505	ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr) ! u'w' at k=0
1506	ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr) ! v'w' at k=0
1507	ts_value(21,sr) = hom(nzb,1,48,sr) ! w"q" at k=0
1508
1509	IF ( ts_value(5,sr) /= 0.0_wp ) THEN
1510	ts_value(22,sr) = ts_value(4,sr)**2 / &
1511	( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
1512	ELSE
1513	ts_value(22,sr) = 10000.0_wp
1514	ENDIF
1515
1516	ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr) ! q*
1517
1518	!
1519	!-- Collect land surface model timeseries
1520	IF ( land_surface ) THEN
1521	ts_value(dots_soil ,sr) = hom(nzb,1,93,sr) ! ghf_eb
1522	ts_value(dots_soil+1,sr) = hom(nzb,1,94,sr) ! shf_eb
1523	ts_value(dots_soil+2,sr) = hom(nzb,1,95,sr) ! qsws_eb
1524	ts_value(dots_soil+3,sr) = hom(nzb,1,96,sr) ! qsws_liq_eb
1525	ts_value(dots_soil+4,sr) = hom(nzb,1,97,sr) ! qsws_soil_eb
1526	ts_value(dots_soil+5,sr) = hom(nzb,1,98,sr) ! qsws_veg_eb
1527	ts_value(dots_soil+6,sr) = hom(nzb,1,99,sr) ! r_a
1528	ts_value(dots_soil+7,sr) = hom(nzb,1,100,sr) ! r_s
1529	ENDIF
1530	!
1531	!-- Collect radiation model timeseries
1532	IF ( radiation ) THEN
1533	ts_value(dots_rad,sr) = hom(nzb,1,101,sr) ! rad_net
1534	ts_value(dots_rad+1,sr) = hom(nzb,1,102,sr) ! rad_lw_in
1535	ts_value(dots_rad+2,sr) = hom(nzb,1,103,sr) ! rad_lw_out
1536	ts_value(dots_rad+3,sr) = hom(nzb,1,104,sr) ! rad_lw_in
1537	ts_value(dots_rad+4,sr) = hom(nzb,1,105,sr) ! rad_lw_out
1538
1539	#if defined ( __rrtmg )
1540	IF ( radiation_scheme == 'rrtmg' ) THEN
1541	ts_value(dots_rad+5,sr) = hom(nzb,1,106,sr) ! rrtm_aldif
1542	ts_value(dots_rad+6,sr) = hom(nzb,1,107,sr) ! rrtm_aldir
1543	ts_value(dots_rad+7,sr) = hom(nzb,1,108,sr) ! rrtm_asdif
1544	ts_value(dots_rad+8,sr) = hom(nzb,1,109,sr) ! rrtm_asdir
1545	ENDIF
1546	#endif
1547
1548	ENDIF
1549
1550	!
1551	!-- Calculate additional statistics provided by the user interface
1552	CALL user_statistics( 'time_series', sr, 0 )
1553
1554	ENDDO ! loop of the subregions
1555
1556	!
1557	!-- If required, sum up horizontal averages for subsequent time averaging
1558	IF ( do_sum ) THEN
1559	IF ( average_count_pr == 0 ) hom_sum = 0.0_wp
1560	hom_sum = hom_sum + hom(:,1,:,:)
1561	average_count_pr = average_count_pr + 1
1562	do_sum = .FALSE.
1563	ENDIF
1564
1565	!
1566	!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
1567	!-- This flag is reset after each time step in time_integration.
1568	flow_statistics_called = .TRUE.
1569
1570	CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
1571
1572
1573	END SUBROUTINE flow_statistics
1574
1575
1576	#else
1577
1578
1579	!------------------------------------------------------------------------------!
1580	! flow statistics - accelerator version
1581	!------------------------------------------------------------------------------!
1582	SUBROUTINE flow_statistics
1583
1584	USE arrays_3d, &
1585	ONLY: ddzu, ddzw, e, hyp, km, kh, nr, p, prho, pt, q, qc, ql, qr, qs, &
1586	qsws, qswst, rho, sa, saswsb, saswst, shf, td_lsa_lpt, td_lsa_q,&
1587	td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us, usws, &
1588	uswst, vsws, v, vg, vpt, vswst, w, w_subs, zw
1589
1590
1591	USE cloud_parameters, &
1592	ONLY: l_d_cp, prr, pt_d_t
1593
1594	USE control_parameters, &
1595	ONLY : average_count_pr, cloud_droplets, cloud_physics, do_sum, &
1596	dt_3d, g, humidity, icloud_scheme, kappa, large_scale_forcing, &
1597	large_scale_subsidence, max_pr_user, message_string, ocean, &
1598	passive_scalar, precipitation, simulated_time, &
1599	use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
1600	ws_scheme_mom, ws_scheme_sca
1601
1602	USE cpulog, &
1603	ONLY: cpu_log, log_point
1604
1605	USE grid_variables, &
1606	ONLY: ddx, ddy
1607
1608	USE indices, &
1609	ONLY: ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums, &
1610	ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzb_diff_s_inner, &
1611	nzb_s_inner, nzt, nzt_diff, rflags_invers
1612
1613	USE kinds
1614
1615	USE land_surface_model_mod, &
1616	ONLY: dots_soil, ghf_eb, land_surface, m_soil, nzb_soil, nzt_soil, &
1617	qsws_eb, qsws_liq_eb, qsws_soil_eb, qsws_veg_eb, r_a, r_s, &
1618	shf_eb, t_soil
1619
1620	USE pegrid
1621
1622	USE radiation_model_mod, &
1623	ONLY: dots_rad, radiation, radiation_scheme, rad_net, &
1624	rad_lw_in, rad_lw_out, rad_sw_in, rad_sw_out
1625
1626	#if defined ( __rrtmg )
1627	USE radiation_model_mod, &
1628	ONLY: rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
1629	#endif
1630
1631	USE statistics
1632
1633	IMPLICIT NONE
1634
1635	INTEGER(iwp) :: i !:
1636	INTEGER(iwp) :: j !:
1637	INTEGER(iwp) :: k !:
1638	INTEGER(iwp) :: nt !:
1639	INTEGER(iwp) :: omp_get_thread_num !:
1640	INTEGER(iwp) :: sr !:
1641	INTEGER(iwp) :: tn !:
1642
1643	LOGICAL :: first !:
1644
1645	REAL(wp) :: dptdz_threshold !:
1646	REAL(wp) :: fac !:
1647	REAL(wp) :: height !:
1648	REAL(wp) :: pts !:
1649	REAL(wp) :: sums_l_eper !:
1650	REAL(wp) :: sums_l_etot !:
1651	REAL(wp) :: s1 !:
1652	REAL(wp) :: s2 !:
1653	REAL(wp) :: s3 !:
1654	REAL(wp) :: s4 !:
1655	REAL(wp) :: s5 !:
1656	REAL(wp) :: s6 !:
1657	REAL(wp) :: s7 !:
1658	REAL(wp) :: ust !:
1659	REAL(wp) :: ust2 !:
1660	REAL(wp) :: u2 !:
1661	REAL(wp) :: vst !:
1662	REAL(wp) :: vst2 !:
1663	REAL(wp) :: v2 !:
1664	REAL(wp) :: w2 !:
1665	REAL(wp) :: z_i(2) !:
1666
1667	REAL(wp) :: dptdz(nzb+1:nzt+1) !:
1668	REAL(wp) :: sums_ll(nzb:nzt+1,2) !:
1669
1670	CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
1671
1672	!
1673	!-- To be on the safe side, check whether flow_statistics has already been
1674	!-- called once after the current time step
1675	IF ( flow_statistics_called ) THEN
1676
1677	message_string = 'flow_statistics is called two times within one ' // &
1678	'timestep'
1679	CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
1680
1681	ENDIF
1682
1683	!$acc data create( sums, sums_l )
1684	!$acc update device( hom )
1685
1686	!
1687	!-- Compute statistics for each (sub-)region
1688	DO sr = 0, statistic_regions
1689
1690	!
1691	!-- Initialize (local) summation array
1692	sums_l = 0.0_wp
1693
1694	!
1695	!-- Store sums that have been computed in other subroutines in summation
1696	!-- array
1697	sums_l(:,11,:) = sums_l_l(:,sr,:) ! mixing length from diffusivities
1698	!-- WARNING: next line still has to be adjusted for OpenMP
1699	sums_l(:,21,0) = sums_wsts_bc_l(:,sr) ! heat flux from advec_s_bc
1700	sums_l(nzb+9,pr_palm,0) = sums_divold_l(sr) ! old divergence from pres
1701	sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr) ! new divergence from pres
1702
1703	!
1704	!-- When calcuating horizontally-averaged total (resolved- plus subgrid-
1705	!-- scale) vertical fluxes and velocity variances by using commonly-
1706	!-- applied Reynolds-based methods ( e.g. <w'pt'> = (w-<w>)*(pt-<pt>) )
1707	!-- in combination with the 5th order advection scheme, pronounced
1708	!-- artificial kinks could be observed in the vertical profiles near the
1709	!-- surface. Please note: these kinks were not related to the model truth,
1710	!-- i.e. these kinks are just related to an evaluation problem.
1711	!-- In order avoid these kinks, vertical fluxes and horizontal as well
1712	!-- vertical velocity variances are calculated directly within the advection
1713	!-- routines, according to the numerical discretization, to evaluate the
1714	!-- statistical quantities as they will appear within the prognostic
1715	!-- equations.
1716	!-- Copy the turbulent quantities, evaluated in the advection routines to
1717	!-- the local array sums_l() for further computations.
1718	IF ( ws_scheme_mom .AND. sr == 0 ) THEN
1719
1720	!
1721	!-- According to the Neumann bc for the horizontal velocity components,
1722	!-- the corresponding fluxes has to satisfiy the same bc.
1723	IF ( ocean ) THEN
1724	sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
1725	sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)
1726	ENDIF
1727
1728	DO i = 0, threads_per_task-1
1729	!
1730	!-- Swap the turbulent quantities evaluated in advec_ws.
1731	sums_l(:,13,i) = sums_wsus_ws_l(:,i) ! wu
1732	sums_l(:,15,i) = sums_wsvs_ws_l(:,i) ! wv
1733	sums_l(:,30,i) = sums_us2_ws_l(:,i) ! u*2
1734	sums_l(:,31,i) = sums_vs2_ws_l(:,i) ! v*2
1735	sums_l(:,32,i) = sums_ws2_ws_l(:,i) ! w*2
1736	sums_l(:,34,i) = sums_l(:,34,i) + 0.5_wp * &
1737	( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) + &
1738	sums_ws2_ws_l(:,i) ) ! e*
1739	DO k = nzb, nzt
1740	sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5_wp * ( &
1741	sums_us2_ws_l(k,i) + &
1742	sums_vs2_ws_l(k,i) + &
1743	sums_ws2_ws_l(k,i) )
1744	ENDDO
1745	ENDDO
1746
1747	ENDIF
1748
1749	IF ( ws_scheme_sca .AND. sr == 0 ) THEN
1750
1751	DO i = 0, threads_per_task-1
1752	sums_l(:,17,i) = sums_wspts_ws_l(:,i) ! wpt from advec_s_ws
1753	IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! wsa
1754	IF ( humidity .OR. passive_scalar ) sums_l(:,49,i) = &
1755	sums_wsqs_ws_l(:,i) !wq
1756	ENDDO
1757
1758	ENDIF
1759	!
1760	!-- Horizontally averaged profiles of horizontal velocities and temperature.
1761	!-- They must have been computed before, because they are already required
1762	!-- for other horizontal averages.
1763	tn = 0
1764
1765	!$OMP PARALLEL PRIVATE( i, j, k, tn )
1766	#if defined( __intel_openmp_bug )
1767	tn = omp_get_thread_num()
1768	#else
1769	!$ tn = omp_get_thread_num()
1770	#endif
1771
1772	!$acc update device( sums_l )
1773
1774	!$OMP DO
1775	!$acc parallel loop gang present( pt, rflags_invers, rmask, sums_l, u, v ) create( s1, s2, s3 )
1776	DO k = nzb, nzt+1
1777	!$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
1778	DO i = nxl, nxr
1779	DO j = nys, nyn
1780	!
1781	!-- k+1 is used in rflags since rflags is set 0 at surface points
1782	s1 = s1 + u(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1783	s2 = s2 + v(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1784	s3 = s3 + pt(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1785	ENDDO
1786	ENDDO
1787	sums_l(k,1,tn) = s1
1788	sums_l(k,2,tn) = s2
1789	sums_l(k,4,tn) = s3
1790	ENDDO
1791	!$acc end parallel loop
1792
1793	!
1794	!-- Horizontally averaged profile of salinity
1795	IF ( ocean ) THEN
1796	!$OMP DO
1797	!$acc parallel loop gang present( rflags_invers, rmask, sums_l, sa ) create( s1 )
1798	DO k = nzb, nzt+1
1799	!$acc loop vector collapse( 2 ) reduction( +: s1 )
1800	DO i = nxl, nxr
1801	DO j = nys, nyn
1802	s1 = s1 + sa(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1803	ENDDO
1804	ENDDO
1805	sums_l(k,23,tn) = s1
1806	ENDDO
1807	!$acc end parallel loop
1808	ENDIF
1809
1810	!
1811	!-- Horizontally averaged profiles of virtual potential temperature,
1812	!-- total water content, specific humidity and liquid water potential
1813	!-- temperature
1814	IF ( humidity ) THEN
1815
1816	!$OMP DO
1817	!$acc parallel loop gang present( q, rflags_invers, rmask, sums_l, vpt ) create( s1, s2 )
1818	DO k = nzb, nzt+1
1819	!$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
1820	DO i = nxl, nxr
1821	DO j = nys, nyn
1822	s1 = s1 + q(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1823	s2 = s2 + vpt(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1824	ENDDO
1825	ENDDO
1826	sums_l(k,41,tn) = s1
1827	sums_l(k,44,tn) = s2
1828	ENDDO
1829	!$acc end parallel loop
1830
1831	IF ( cloud_physics ) THEN
1832	!$OMP DO
1833	!$acc parallel loop gang present( pt, q, ql, rflags_invers, rmask, sums_l ) create( s1, s2 )
1834	DO k = nzb, nzt+1
1835	!$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
1836	DO i = nxl, nxr
1837	DO j = nys, nyn
1838	s1 = s1 + ( q(k,j,i) - ql(k,j,i) ) * &
1839	rmask(j,i,sr) * rflags_invers(j,i,k+1)
1840	s2 = s2 + ( pt(k,j,i) + l_d_cppt_d_t(k) ql(k,j,i) ) * &
1841	rmask(j,i,sr) * rflags_invers(j,i,k+1)
1842	ENDDO
1843	ENDDO
1844	sums_l(k,42,tn) = s1
1845	sums_l(k,43,tn) = s2
1846	ENDDO
1847	!$acc end parallel loop
1848	ENDIF
1849	ENDIF
1850
1851	!
1852	!-- Horizontally averaged profiles of passive scalar
1853	IF ( passive_scalar ) THEN
1854	!$OMP DO
1855	!$acc parallel loop gang present( q, rflags_invers, rmask, sums_l ) create( s1 )
1856	DO k = nzb, nzt+1
1857	!$acc loop vector collapse( 2 ) reduction( +: s1 )
1858	DO i = nxl, nxr
1859	DO j = nys, nyn
1860	s1 = s1 + q(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1861	ENDDO
1862	ENDDO
1863	sums_l(k,41,tn) = s1
1864	ENDDO
1865	!$acc end parallel loop
1866	ENDIF
1867	!$OMP END PARALLEL
1868
1869	!
1870	!-- Summation of thread sums
1871	IF ( threads_per_task > 1 ) THEN
1872	DO i = 1, threads_per_task-1
1873	!$acc parallel present( sums_l )
1874	sums_l(:,1,0) = sums_l(:,1,0) + sums_l(:,1,i)
1875	sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i)
1876	sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i)
1877	!$acc end parallel
1878	IF ( ocean ) THEN
1879	!$acc parallel present( sums_l )
1880	sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)
1881	!$acc end parallel
1882	ENDIF
1883	IF ( humidity ) THEN
1884	!$acc parallel present( sums_l )
1885	sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
1886	sums_l(:,44,0) = sums_l(:,44,0) + sums_l(:,44,i)
1887	!$acc end parallel
1888	IF ( cloud_physics ) THEN
1889	!$acc parallel present( sums_l )
1890	sums_l(:,42,0) = sums_l(:,42,0) + sums_l(:,42,i)
1891	sums_l(:,43,0) = sums_l(:,43,0) + sums_l(:,43,i)
1892	!$acc end parallel
1893	ENDIF
1894	ENDIF
1895	IF ( passive_scalar ) THEN
1896	!$acc parallel present( sums_l )
1897	sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
1898	!$acc end parallel
1899	ENDIF
1900	ENDDO
1901	ENDIF
1902
1903	#if defined( __parallel )
1904	!
1905	!-- Compute total sum from local sums
1906	!$acc update host( sums_l )
1907	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1908	CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL, &
1909	MPI_SUM, comm2d, ierr )
1910	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1911	CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL, &
1912	MPI_SUM, comm2d, ierr )
1913	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1914	CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL, &
1915	MPI_SUM, comm2d, ierr )
1916	IF ( ocean ) THEN
1917	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1918	CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb, &
1919	MPI_REAL, MPI_SUM, comm2d, ierr )
1920	ENDIF
1921	IF ( humidity ) THEN
1922	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1923	CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb, &
1924	MPI_REAL, MPI_SUM, comm2d, ierr )
1925	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1926	CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
1927	MPI_REAL, MPI_SUM, comm2d, ierr )
1928	IF ( cloud_physics ) THEN
1929	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1930	CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb, &
1931	MPI_REAL, MPI_SUM, comm2d, ierr )
1932	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1933	CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb, &
1934	MPI_REAL, MPI_SUM, comm2d, ierr )
1935	ENDIF
1936	ENDIF
1937
1938	IF ( passive_scalar ) THEN
1939	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1940	CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
1941	MPI_REAL, MPI_SUM, comm2d, ierr )
1942	ENDIF
1943	!$acc update device( sums )
1944	#else
1945	!$acc parallel present( sums, sums_l )
1946	sums(:,1) = sums_l(:,1,0)
1947	sums(:,2) = sums_l(:,2,0)
1948	sums(:,4) = sums_l(:,4,0)
1949	!$acc end parallel
1950	IF ( ocean ) THEN
1951	!$acc parallel present( sums, sums_l )
1952	sums(:,23) = sums_l(:,23,0)
1953	!$acc end parallel
1954	ENDIF
1955	IF ( humidity ) THEN
1956	!$acc parallel present( sums, sums_l )
1957	sums(:,44) = sums_l(:,44,0)
1958	sums(:,41) = sums_l(:,41,0)
1959	!$acc end parallel
1960	IF ( cloud_physics ) THEN
1961	!$acc parallel present( sums, sums_l )
1962	sums(:,42) = sums_l(:,42,0)
1963	sums(:,43) = sums_l(:,43,0)
1964	!$acc end parallel
1965	ENDIF
1966	ENDIF
1967	IF ( passive_scalar ) THEN
1968	!$acc parallel present( sums, sums_l )
1969	sums(:,41) = sums_l(:,41,0)
1970	!$acc end parallel
1971	ENDIF
1972	#endif
1973
1974	!
1975	!-- Final values are obtained by division by the total number of grid points
1976	!-- used for summation. After that store profiles.
1977	!$acc parallel present( hom, ngp_2dh, ngp_2dh_s_inner, sums )
1978	sums(:,1) = sums(:,1) / ngp_2dh(sr)
1979	sums(:,2) = sums(:,2) / ngp_2dh(sr)
1980	sums(:,4) = sums(:,4) / ngp_2dh_s_inner(:,sr)
1981	hom(:,1,1,sr) = sums(:,1) ! u
1982	hom(:,1,2,sr) = sums(:,2) ! v
1983	hom(:,1,4,sr) = sums(:,4) ! pt
1984	!$acc end parallel
1985
1986	!
1987	!-- Salinity
1988	IF ( ocean ) THEN
1989	!$acc parallel present( hom, ngp_2dh_s_inner, sums )
1990	sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr)
1991	hom(:,1,23,sr) = sums(:,23) ! sa
1992	!$acc end parallel
1993	ENDIF
1994
1995	!
1996	!-- Humidity and cloud parameters
1997	IF ( humidity ) THEN
1998	!$acc parallel present( hom, ngp_2dh_s_inner, sums )
1999	sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr)
2000	sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
2001	hom(:,1,44,sr) = sums(:,44) ! vpt
2002	hom(:,1,41,sr) = sums(:,41) ! qv (q)
2003	!$acc end parallel
2004	IF ( cloud_physics ) THEN
2005	!$acc parallel present( hom, ngp_2dh_s_inner, sums )
2006	sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr)
2007	sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr)
2008	hom(:,1,42,sr) = sums(:,42) ! qv
2009	hom(:,1,43,sr) = sums(:,43) ! pt
2010	!$acc end parallel
2011	ENDIF
2012	ENDIF
2013
2014	!
2015	!-- Passive scalar
2016	IF ( passive_scalar ) THEN
2017	!$acc parallel present( hom, ngp_2dh_s_inner, sums )
2018	sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
2019	hom(:,1,41,sr) = sums(:,41) ! s (q)
2020	!$acc end parallel
2021	ENDIF
2022
2023	!
2024	!-- Horizontally averaged profiles of the remaining prognostic variables,
2025	!-- variances, the total and the perturbation energy (single values in last
2026	!-- column of sums_l) and some diagnostic quantities.
2027	!-- NOTE: for simplicity, nzb_s_inner is used below, although strictly
2028	!-- ---- speaking the following k-loop would have to be split up and
2029	!-- rearranged according to the staggered grid.
2030	!-- However, this implies no error since staggered velocity components
2031	!-- are zero at the walls and inside buildings.
2032	tn = 0
2033	#if defined( __intel_openmp_bug )
2034	!$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, &
2035	!$OMP tn, ust, ust2, u2, vst, vst2, v2, w2 )
2036	tn = omp_get_thread_num()
2037	#else
2038	!$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2, w2 )
2039	!$ tn = omp_get_thread_num()
2040	#endif
2041	!$OMP DO
2042	!$acc parallel loop gang present( e, hom, kh, km, p, pt, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3, s4, s5, s6, s7 )
2043	DO k = nzb, nzt+1
2044	!$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5, s6, s7 )
2045	DO i = nxl, nxr
2046	DO j = nys, nyn
2047	!
2048	!-- Prognostic and diagnostic variables
2049	s1 = s1 + w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2050	s2 = s2 + e(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2051	s3 = s3 + km(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2052	s4 = s4 + kh(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2053	s5 = s5 + p(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2054	s6 = s6 + ( pt(k,j,i)-hom(k,1,4,sr) )*2 rmask(j,i,sr) * &
2055	rflags_invers(j,i,k+1)
2056	!
2057	!-- Higher moments
2058	!-- (Computation of the skewness of w further below)
2059	s7 = s7 + w(k,j,i)*3 rmask(j,i,sr) * rflags_invers(j,i,k+1)
2060	ENDDO
2061	ENDDO
2062	sums_l(k,3,tn) = s1
2063	sums_l(k,8,tn) = s2
2064	sums_l(k,9,tn) = s3
2065	sums_l(k,10,tn) = s4
2066	sums_l(k,40,tn) = s5
2067	sums_l(k,33,tn) = s6
2068	sums_l(k,38,tn) = s7
2069	ENDDO
2070	!$acc end parallel loop
2071
2072	IF ( humidity ) THEN
2073	!$OMP DO
2074	!$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l ) create( s1 )
2075	DO k = nzb, nzt+1
2076	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2077	DO i = nxl, nxr
2078	DO j = nys, nyn
2079	s1 = s1 + ( q(k,j,i)-hom(k,1,41,sr) )*2 rmask(j,i,sr) * &
2080	rflags_invers(j,i,k+1)
2081	ENDDO
2082	ENDDO
2083	sums_l(k,70,tn) = s1
2084	ENDDO
2085	!$acc end parallel loop
2086	ENDIF
2087
2088	!
2089	!-- Total and perturbation energy for the total domain (being
2090	!-- collected in the last column of sums_l).
2091	!$OMP DO
2092	!$acc parallel loop collapse(3) present( rflags_invers, rmask, u, v, w ) reduction(+:s1)
2093	DO i = nxl, nxr
2094	DO j = nys, nyn
2095	DO k = nzb, nzt+1
2096	s1 = s1 + 0.5_wp * &
2097	( u(k,j,i)2 + v(k,j,i)2 + w(k,j,i)*2 ) &
2098	rmask(j,i,sr) * rflags_invers(j,i,k+1)
2099	ENDDO
2100	ENDDO
2101	ENDDO
2102	!$acc end parallel loop
2103	!$acc parallel present( sums_l )
2104	sums_l(nzb+4,pr_palm,tn) = s1
2105	!$acc end parallel
2106
2107	!$OMP DO
2108	!$acc parallel present( rmask, sums_l, us, usws, vsws, ts ) create( s1, s2, s3, s4 )
2109	!$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4 )
2110	DO i = nxl, nxr
2111	DO j = nys, nyn
2112	!
2113	!-- 2D-arrays (being collected in the last column of sums_l)
2114	s1 = s1 + us(j,i) * rmask(j,i,sr)
2115	s2 = s2 + usws(j,i) * rmask(j,i,sr)
2116	s3 = s3 + vsws(j,i) * rmask(j,i,sr)
2117	s4 = s4 + ts(j,i) * rmask(j,i,sr)
2118	ENDDO
2119	ENDDO
2120	sums_l(nzb,pr_palm,tn) = s1
2121	sums_l(nzb+1,pr_palm,tn) = s2
2122	sums_l(nzb+2,pr_palm,tn) = s3
2123	sums_l(nzb+3,pr_palm,tn) = s4
2124	!$acc end parallel
2125
2126	IF ( humidity ) THEN
2127	!$acc parallel present( qs, rmask, sums_l ) create( s1 )
2128	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2129	DO i = nxl, nxr
2130	DO j = nys, nyn
2131	s1 = s1 + qs(j,i) * rmask(j,i,sr)
2132	ENDDO
2133	ENDDO
2134	sums_l(nzb+12,pr_palm,tn) = s1
2135	!$acc end parallel
2136	ENDIF
2137
2138	!
2139	!-- Computation of statistics when ws-scheme is not used. Else these
2140	!-- quantities are evaluated in the advection routines.
2141	IF ( .NOT. ws_scheme_mom .OR. sr /= 0 ) THEN
2142
2143	!$OMP DO
2144	!$acc parallel loop gang present( u, v, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3, s4, ust2, vst2, w2 )
2145	DO k = nzb, nzt+1
2146	!$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4 )
2147	DO i = nxl, nxr
2148	DO j = nys, nyn
2149	ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
2150	vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
2151	w2 = w(k,j,i)**2
2152
2153	s1 = s1 + ust2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2154	s2 = s2 + vst2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2155	s3 = s3 + w2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2156	!
2157	!-- Perturbation energy
2158	s4 = s4 + 0.5_wp * ( ust2 + vst2 + w2 ) * rmask(j,i,sr) * &
2159	rflags_invers(j,i,k+1)
2160	ENDDO
2161	ENDDO
2162	sums_l(k,30,tn) = s1
2163	sums_l(k,31,tn) = s2
2164	sums_l(k,32,tn) = s3
2165	sums_l(k,34,tn) = s4
2166	ENDDO
2167	!$acc end parallel loop
2168	!
2169	!-- Total perturbation TKE
2170	!$OMP DO
2171	!$acc parallel present( sums_l ) create( s1 )
2172	!$acc loop reduction( +: s1 )
2173	DO k = nzb, nzt+1
2174	s1 = s1 + sums_l(k,34,tn)
2175	ENDDO
2176	sums_l(nzb+5,pr_palm,tn) = s1
2177	!$acc end parallel
2178
2179	ENDIF
2180
2181	!
2182	!-- Horizontally averaged profiles of the vertical fluxes
2183
2184	!
2185	!-- Subgridscale fluxes.
2186	!-- WARNING: If a Prandtl-layer is used (k=nzb for flat terrain), the fluxes
2187	!-- ------- should be calculated there in a different way. This is done
2188	!-- in the next loop further below, where results from this loop are
2189	!-- overwritten. However, THIS WORKS IN CASE OF FLAT TERRAIN ONLY!
2190	!-- The non-flat case still has to be handled.
2191	!-- NOTE: for simplicity, nzb_s_inner is used below, although
2192	!-- ---- strictly speaking the following k-loop would have to be
2193	!-- split up according to the staggered grid.
2194	!-- However, this implies no error since staggered velocity
2195	!-- components are zero at the walls and inside buildings.
2196	!$OMP DO
2197	!$acc parallel loop gang present( ddzu, kh, km, pt, u, v, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3 )
2198	DO k = nzb, nzt_diff
2199	!$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
2200	DO i = nxl, nxr
2201	DO j = nys, nyn
2202
2203	!
2204	!-- Momentum flux w"u"
2205	s1 = s1 - 0.25_wp * ( &
2206	km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
2207	) * ( &
2208	( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
2209	+ ( w(k,j,i) - w(k,j,i-1) ) * ddx &
2210	) &
2211	* rmask(j,i,sr) * rflags_invers(j,i,k+1)
2212	!
2213	!-- Momentum flux w"v"
2214	s2 = s2 - 0.25_wp * ( &
2215	km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
2216	) * ( &
2217	( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
2218	+ ( w(k,j,i) - w(k,j-1,i) ) * ddy &
2219	) &
2220	* rmask(j,i,sr) * rflags_invers(j,i,k+1)
2221	!
2222	!-- Heat flux w"pt"
2223	s3 = s3 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) ) &
2224	* ( pt(k+1,j,i) - pt(k,j,i) ) &
2225	* ddzu(k+1) * rmask(j,i,sr) &
2226	* rflags_invers(j,i,k+1)
2227	ENDDO
2228	ENDDO
2229	sums_l(k,12,tn) = s1
2230	sums_l(k,14,tn) = s2
2231	sums_l(k,16,tn) = s3
2232	ENDDO
2233	!$acc end parallel loop
2234
2235	!
2236	!-- Salinity flux w"sa"
2237	IF ( ocean ) THEN
2238	!$acc parallel loop gang present( ddzu, kh, sa, rflags_invers, rmask, sums_l ) create( s1 )
2239	DO k = nzb, nzt_diff
2240	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2241	DO i = nxl, nxr
2242	DO j = nys, nyn
2243	s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) ) &
2244	* ( sa(k+1,j,i) - sa(k,j,i) ) &
2245	* ddzu(k+1) * rmask(j,i,sr) &
2246	* rflags_invers(j,i,k+1)
2247	ENDDO
2248	ENDDO
2249	sums_l(k,65,tn) = s1
2250	ENDDO
2251	!$acc end parallel loop
2252	ENDIF
2253
2254	!
2255	!-- Buoyancy flux, water flux (humidity flux) w"q"
2256	IF ( humidity ) THEN
2257
2258	!$acc parallel loop gang present( ddzu, kh, q, vpt, rflags_invers, rmask, sums_l ) create( s1, s2 )
2259	DO k = nzb, nzt_diff
2260	!$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2261	DO i = nxl, nxr
2262	DO j = nys, nyn
2263	s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) ) &
2264	* ( vpt(k+1,j,i) - vpt(k,j,i) ) &
2265	* ddzu(k+1) * rmask(j,i,sr) &
2266	* rflags_invers(j,i,k+1)
2267	s2 = s2 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) ) &
2268	* ( q(k+1,j,i) - q(k,j,i) ) &
2269	* ddzu(k+1) * rmask(j,i,sr) &
2270	* rflags_invers(j,i,k+1)
2271	ENDDO
2272	ENDDO
2273	sums_l(k,45,tn) = s1
2274	sums_l(k,48,tn) = s2
2275	ENDDO
2276	!$acc end parallel loop
2277
2278	IF ( cloud_physics ) THEN
2279
2280	!$acc parallel loop gang present( ddzu, kh, q, ql, rflags_invers, rmask, sums_l ) create( s1 )
2281	DO k = nzb, nzt_diff
2282	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2283	DO i = nxl, nxr
2284	DO j = nys, nyn
2285	s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) ) &
2286	* ( ( q(k+1,j,i) - ql(k+1,j,i) ) &
2287	- ( q(k,j,i) - ql(k,j,i) ) ) &
2288	* ddzu(k+1) * rmask(j,i,sr) &
2289	* rflags_invers(j,i,k+1)
2290	ENDDO
2291	ENDDO
2292	sums_l(k,51,tn) = s1
2293	ENDDO
2294	!$acc end parallel loop
2295
2296	ENDIF
2297
2298	ENDIF
2299	!
2300	!-- Passive scalar flux
2301	IF ( passive_scalar ) THEN
2302
2303	!$acc parallel loop gang present( ddzu, kh, q, rflags_invers, rmask, sums_l ) create( s1 )
2304	DO k = nzb, nzt_diff
2305	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2306	DO i = nxl, nxr
2307	DO j = nys, nyn
2308	s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) ) &
2309	* ( q(k+1,j,i) - q(k,j,i) ) &
2310	* ddzu(k+1) * rmask(j,i,sr) &
2311	* rflags_invers(j,i,k+1)
2312	ENDDO
2313	ENDDO
2314	sums_l(k,48,tn) = s1
2315	ENDDO
2316	!$acc end parallel loop
2317
2318	ENDIF
2319
2320	IF ( use_surface_fluxes ) THEN
2321
2322	!$OMP DO
2323	!$acc parallel present( rmask, shf, sums_l, usws, vsws ) create( s1, s2, s3, s4, s5 )
2324	!$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5 )
2325	DO i = nxl, nxr
2326	DO j = nys, nyn
2327	!
2328	!-- Subgridscale fluxes in the Prandtl layer
2329	s1 = s1 + usws(j,i) * rmask(j,i,sr) ! w"u"
2330	s2 = s2 + vsws(j,i) * rmask(j,i,sr) ! w"v"
2331	s3 = s3 + shf(j,i) * rmask(j,i,sr) ! w"pt"
2332	s4 = s4 + 0.0_wp * rmask(j,i,sr) ! u"pt"
2333	s5 = s5 + 0.0_wp * rmask(j,i,sr) ! v"pt"
2334	ENDDO
2335	ENDDO
2336	sums_l(nzb,12,tn) = s1
2337	sums_l(nzb,14,tn) = s2
2338	sums_l(nzb,16,tn) = s3
2339	sums_l(nzb,58,tn) = s4
2340	sums_l(nzb,61,tn) = s5
2341	!$acc end parallel
2342
2343	IF ( ocean ) THEN
2344
2345	!$OMP DO
2346	!$acc parallel present( rmask, saswsb, sums_l ) create( s1 )
2347	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2348	DO i = nxl, nxr
2349	DO j = nys, nyn
2350	s1 = s1 + saswsb(j,i) * rmask(j,i,sr) ! w"sa"
2351	ENDDO
2352	ENDDO
2353	sums_l(nzb,65,tn) = s1
2354	!$acc end parallel
2355
2356	ENDIF
2357
2358	IF ( humidity ) THEN
2359
2360	!$OMP DO
2361	!$acc parallel present( pt, q, qsws, rmask, shf, sums_l ) create( s1, s2 )
2362	!$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2363	DO i = nxl, nxr
2364	DO j = nys, nyn
2365	s1 = s1 + qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
2366	s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) * shf(j,i) &
2367	+ 0.61_wp * pt(nzb,j,i) * qsws(j,i) )
2368	ENDDO
2369	ENDDO
2370	sums_l(nzb,48,tn) = s1
2371	sums_l(nzb,45,tn) = s2
2372	!$acc end parallel
2373
2374	IF ( cloud_droplets ) THEN
2375
2376	!$OMP DO
2377	!$acc parallel present( pt, q, ql, qsws, rmask, shf, sums_l ) create( s1 )
2378	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2379	DO i = nxl, nxr
2380	DO j = nys, nyn
2381	s1 = s1 + ( ( 1.0_wp + &
2382	0.61_wp * q(nzb,j,i) - ql(nzb,j,i) ) * &
2383	shf(j,i) + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )
2384	ENDDO
2385	ENDDO
2386	sums_l(nzb,45,tn) = s1
2387	!$acc end parallel
2388
2389	ENDIF
2390
2391	IF ( cloud_physics ) THEN
2392
2393	!$OMP DO
2394	!$acc parallel present( qsws, rmask, sums_l ) create( s1 )
2395	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2396	DO i = nxl, nxr
2397	DO j = nys, nyn
2398	!
2399	!-- Formula does not work if ql(nzb) /= 0.0
2400	s1 = s1 + qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
2401	ENDDO
2402	ENDDO
2403	sums_l(nzb,51,tn) = s1
2404	!$acc end parallel
2405
2406	ENDIF
2407
2408	ENDIF
2409
2410	IF ( passive_scalar ) THEN
2411
2412	!$OMP DO
2413	!$acc parallel present( qsws, rmask, sums_l ) create( s1 )
2414	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2415	DO i = nxl, nxr
2416	DO j = nys, nyn
2417	s1 = s1 + qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
2418	ENDDO
2419	ENDDO
2420	sums_l(nzb,48,tn) = s1
2421	!$acc end parallel
2422
2423	ENDIF
2424
2425	ENDIF
2426
2427	!
2428	!-- Subgridscale fluxes at the top surface
2429	IF ( use_top_fluxes ) THEN
2430
2431	!$OMP DO
2432	!$acc parallel present( rmask, sums_l, tswst, uswst, vswst ) create( s1, s2, s3, s4, s5 )
2433	!$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5 )
2434	DO i = nxl, nxr
2435	DO j = nys, nyn
2436	s1 = s1 + uswst(j,i) * rmask(j,i,sr) ! w"u"
2437	s2 = s2 + vswst(j,i) * rmask(j,i,sr) ! w"v"
2438	s3 = s3 + tswst(j,i) * rmask(j,i,sr) ! w"pt"
2439	s4 = s4 + 0.0_wp * rmask(j,i,sr) ! u"pt"
2440	s5 = s5 + 0.0_wp * rmask(j,i,sr) ! v"pt"
2441	ENDDO
2442	ENDDO
2443	sums_l(nzt:nzt+1,12,tn) = s1
2444	sums_l(nzt:nzt+1,14,tn) = s2
2445	sums_l(nzt:nzt+1,16,tn) = s3
2446	sums_l(nzt:nzt+1,58,tn) = s4
2447	sums_l(nzt:nzt+1,61,tn) = s5
2448	!$acc end parallel
2449
2450	IF ( ocean ) THEN
2451
2452	!$OMP DO
2453	!$acc parallel present( rmask, saswst, sums_l ) create( s1 )
2454	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2455	DO i = nxl, nxr
2456	DO j = nys, nyn
2457	s1 = s1 + saswst(j,i) * rmask(j,i,sr) ! w"sa"
2458	ENDDO
2459	ENDDO
2460	sums_l(nzt,65,tn) = s1
2461	!$acc end parallel
2462
2463	ENDIF
2464
2465	IF ( humidity ) THEN
2466
2467	!$OMP DO
2468	!$acc parallel present( pt, q, qswst, rmask, tswst, sums_l ) create( s1, s2 )
2469	!$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2470	DO i = nxl, nxr
2471	DO j = nys, nyn
2472	s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
2473	s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) * tswst(j,i) +&
2474	0.61_wp * pt(nzt,j,i) * qswst(j,i) )
2475	ENDDO
2476	ENDDO
2477	sums_l(nzt,48,tn) = s1
2478	sums_l(nzt,45,tn) = s2
2479	!$acc end parallel
2480
2481	IF ( cloud_droplets ) THEN
2482
2483	!$OMP DO
2484	!$acc parallel present( pt, q, ql, qswst, rmask, tswst, sums_l ) create( s1 )
2485	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2486	DO i = nxl, nxr
2487	DO j = nys, nyn
2488	s1 = s1 + ( ( 1.0_wp + &
2489	0.61_wp * q(nzt,j,i) - ql(nzt,j,i) ) * &
2490	tswst(j,i) + &
2491	0.61_wp * pt(nzt,j,i) * qswst(j,i) )
2492	ENDDO
2493	ENDDO
2494	sums_l(nzt,45,tn) = s1
2495	!$acc end parallel
2496
2497	ENDIF
2498
2499	IF ( cloud_physics ) THEN
2500
2501	!$OMP DO
2502	!$acc parallel present( qswst, rmask, sums_l ) create( s1 )
2503	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2504	DO i = nxl, nxr
2505	DO j = nys, nyn
2506	!
2507	!-- Formula does not work if ql(nzb) /= 0.0
2508	s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
2509	ENDDO
2510	ENDDO
2511	sums_l(nzt,51,tn) = s1
2512	!$acc end parallel
2513
2514	ENDIF
2515
2516	ENDIF
2517
2518	IF ( passive_scalar ) THEN
2519
2520	!$OMP DO
2521	!$acc parallel present( qswst, rmask, sums_l ) create( s1 )
2522	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2523	DO i = nxl, nxr
2524	DO j = nys, nyn
2525	s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
2526	ENDDO
2527	ENDDO
2528	sums_l(nzt,48,tn) = s1
2529	!$acc end parallel
2530
2531	ENDIF
2532
2533	ENDIF
2534
2535	!
2536	!-- Resolved fluxes (can be computed for all horizontal points)
2537	!-- NOTE: for simplicity, nzb_s_inner is used below, although strictly
2538	!-- ---- speaking the following k-loop would have to be split up and
2539	!-- rearranged according to the staggered grid.
2540	!$acc parallel loop gang present( hom, pt, rflags_invers, rmask, sums_l, u, v, w ) create( s1, s2, s3 )
2541	DO k = nzb, nzt_diff
2542	!$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
2543	DO i = nxl, nxr
2544	DO j = nys, nyn
2545	ust = 0.5_wp * ( u(k,j,i) - hom(k,1,1,sr) + &
2546	u(k+1,j,i) - hom(k+1,1,1,sr) )
2547	vst = 0.5_wp * ( v(k,j,i) - hom(k,1,2,sr) + &
2548	v(k+1,j,i) - hom(k+1,1,2,sr) )
2549	pts = 0.5_wp * ( pt(k,j,i) - hom(k,1,4,sr) + &
2550	pt(k+1,j,i) - hom(k+1,1,4,sr) )
2551	!
2552	!-- Higher moments
2553	s1 = s1 + pts * w(k,j,i)*2 rmask(j,i,sr) * rflags_invers(j,i,k+1)
2554	s2 = s2 + pts*2 w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2555	!
2556	!-- Energy flux we (has to be adjusted?)
2557	s3 = s3 + w(k,j,i) * 0.5_wp * ( ust2 + vst2 + w(k,j,i)**2 )&
2558	* rmask(j,i,sr) * rflags_invers(j,i,k+1)
2559	ENDDO
2560	ENDDO
2561	sums_l(k,35,tn) = s1
2562	sums_l(k,36,tn) = s2
2563	sums_l(k,37,tn) = s3
2564	ENDDO
2565	!$acc end parallel loop
2566
2567	!
2568	!-- Salinity flux and density (density does not belong to here,
2569	!-- but so far there is no other suitable place to calculate)
2570	IF ( ocean ) THEN
2571
2572	IF( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN
2573
2574	!$acc parallel loop gang present( hom, rflags_invers, rmask, sa, sums_l, w ) create( s1 )
2575	DO k = nzb, nzt_diff
2576	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2577	DO i = nxl, nxr
2578	DO j = nys, nyn
2579	s1 = s1 + 0.5_wp * ( sa(k,j,i) - hom(k,1,23,sr) + &
2580	sa(k+1,j,i) - hom(k+1,1,23,sr) ) &
2581	* w(k,j,i) * rmask(j,i,sr) &
2582	* rflags_invers(j,i,k+1)
2583	ENDDO
2584	ENDDO
2585	sums_l(k,66,tn) = s1
2586	ENDDO
2587	!$acc end parallel loop
2588
2589	ENDIF
2590
2591	!$acc parallel loop gang present( rflags_invers, rho, prho, rmask, sums_l ) create( s1, s2 )
2592	DO k = nzb, nzt_diff
2593	!$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2594	DO i = nxl, nxr
2595	DO j = nys, nyn
2596	s1 = s1 + rho(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2597	s2 = s2 + prho(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2598	ENDDO
2599	ENDDO
2600	sums_l(k,64,tn) = s1
2601	sums_l(k,71,tn) = s2
2602	ENDDO
2603	!$acc end parallel loop
2604
2605	ENDIF
2606
2607	!
2608	!-- Buoyancy flux, water flux, humidity flux, liquid water
2609	!-- content, rain drop concentration and rain water content
2610	IF ( humidity ) THEN
2611
2612	IF ( cloud_physics .OR. cloud_droplets ) THEN
2613
2614	!$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, vpt, w ) create( s1 )
2615	DO k = nzb, nzt_diff
2616	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2617	DO i = nxl, nxr
2618	DO j = nys, nyn
2619	s1 = s1 + 0.5_wp * ( vpt(k,j,i) - hom(k,1,44,sr) + &
2620	vpt(k+1,j,i) - hom(k+1,1,44,sr) ) * &
2621	w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2622	ENDDO
2623	ENDDO
2624	sums_l(k,46,tn) = s1
2625	ENDDO
2626	!$acc end parallel loop
2627
2628	IF ( .NOT. cloud_droplets ) THEN
2629
2630	!$acc parallel loop gang present( hom, q, ql, rflags_invers, rmask, sums_l, w ) create( s1 )
2631	DO k = nzb, nzt_diff
2632	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2633	DO i = nxl, nxr
2634	DO j = nys, nyn
2635	s1 = s1 + 0.5_wp * ( ( q(k,j,i) - ql(k,j,i) ) - hom(k,1,42,sr) + &
2636	( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) ) &
2637	* w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2638	ENDDO
2639	ENDDO
2640	sums_l(k,52,tn) = s1
2641	ENDDO
2642	!$acc end parallel loop
2643
2644	IF ( icloud_scheme == 0 ) THEN
2645
2646	!$acc parallel loop gang present( qc, ql, rflags_invers, rmask, sums_l ) create( s1, s2 )
2647	DO k = nzb, nzt_diff
2648	!$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2649	DO i = nxl, nxr
2650	DO j = nys, nyn
2651	s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2652	s2 = s2 + qc(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2653	ENDDO
2654	ENDDO
2655	sums_l(k,54,tn) = s1
2656	sums_l(k,75,tn) = s2
2657	ENDDO
2658	!$acc end parallel loop
2659
2660	IF ( precipitation ) THEN
2661
2662	!$acc parallel loop gang present( nr, qr, prr, rflags_invers, rmask, sums_l ) create( s1, s2, s3 )
2663	DO k = nzb, nzt_diff
2664	!$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
2665	DO i = nxl, nxr
2666	DO j = nys, nyn
2667	s1 = s1 + nr(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2668	s2 = s2 + qr(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2669	s3 = s3 + prr(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2670	ENDDO
2671	ENDDO
2672	sums_l(k,73,tn) = s1
2673	sums_l(k,74,tn) = s2
2674	sums_l(k,76,tn) = s3
2675	ENDDO
2676	!$acc end parallel loop
2677
2678	ENDIF
2679
2680	ELSE
2681
2682	!$acc parallel loop gang present( ql, rflags_invers, rmask, sums_l ) create( s1 )
2683	DO k = nzb, nzt_diff
2684	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2685	DO i = nxl, nxr
2686	DO j = nys, nyn
2687	s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2688	ENDDO
2689	ENDDO
2690	sums_l(k,54,tn) = s1
2691	ENDDO
2692	!$acc end parallel loop
2693
2694	ENDIF
2695
2696	ELSE
2697
2698	!$acc parallel loop gang present( ql, rflags_invers, rmask, sums_l ) create( s1 )
2699	DO k = nzb, nzt_diff
2700	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2701	DO i = nxl, nxr
2702	DO j = nys, nyn
2703	s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2704	ENDDO
2705	ENDDO
2706	sums_l(k,54,tn) = s1
2707	ENDDO
2708	!$acc end parallel loop
2709
2710	ENDIF
2711
2712	ELSE
2713
2714	IF( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN
2715
2716	!$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, vpt, w ) create( s1 )
2717	DO k = nzb, nzt_diff
2718	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2719	DO i = nxl, nxr
2720	DO j = nys, nyn
2721	s1 = s1 + 0.5_wp * ( vpt(k,j,i) - hom(k,1,44,sr) + &
2722	vpt(k+1,j,i) - hom(k+1,1,44,sr) ) &
2723	* w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2724	ENDDO
2725	ENDDO
2726	sums_l(k,46,tn) = s1
2727	ENDDO
2728	!$acc end parallel loop
2729
2730	ELSEIF ( ws_scheme_sca .AND. sr == 0 ) THEN
2731
2732	!$acc parallel loop present( hom, sums_l )
2733	DO k = nzb, nzt_diff
2734	sums_l(k,46,tn) = ( 1.0_wp + 0.61_wp * hom(k,1,41,sr) ) * sums_l(k,17,tn) + &
2735	0.61_wp * hom(k,1,4,sr) * sums_l(k,49,tn)
2736	ENDDO
2737	!$acc end parallel loop
2738
2739	ENDIF
2740
2741	ENDIF
2742
2743	ENDIF
2744	!
2745	!-- Passive scalar flux
2746	IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca .OR. sr /= 0 ) ) THEN
2747
2748	!$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l, w ) create( s1 )
2749	DO k = nzb, nzt_diff
2750	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2751	DO i = nxl, nxr
2752	DO j = nys, nyn
2753	s1 = s1 + 0.5_wp * ( q(k,j,i) - hom(k,1,41,sr) + &
2754	q(k+1,j,i) - hom(k+1,1,41,sr) ) &
2755	* w(k,j,i) * rmask(j,i,sr) &
2756	* rflags_invers(j,i,k+1)
2757	ENDDO
2758	ENDDO
2759	sums_l(k,49,tn) = s1
2760	ENDDO
2761	!$acc end parallel loop
2762
2763	ENDIF
2764
2765	!
2766	!-- For speed optimization fluxes which have been computed in part directly
2767	!-- inside the WS advection routines are treated seperatly
2768	!-- Momentum fluxes first:
2769	IF ( .NOT. ws_scheme_mom .OR. sr /= 0 ) THEN
2770
2771	!$OMP DO
2772	!$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, u, v, w ) create( s1, s2 )
2773	DO k = nzb, nzt_diff
2774	!$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2775	DO i = nxl, nxr
2776	DO j = nys, nyn
2777	ust = 0.5_wp * ( u(k,j,i) - hom(k,1,1,sr) + &
2778	u(k+1,j,i) - hom(k+1,1,1,sr) )
2779	vst = 0.5_wp * ( v(k,j,i) - hom(k,1,2,sr) + &
2780	v(k+1,j,i) - hom(k+1,1,2,sr) )
2781	!
2782	!-- Momentum flux wu
2783	s1 = s1 + 0.5_wp * ( w(k,j,i-1) + w(k,j,i) ) &
2784	* ust * rmask(j,i,sr) &
2785	* rflags_invers(j,i,k+1)
2786	!
2787	!-- Momentum flux wv
2788	s2 = s2 + 0.5_wp * ( w(k,j-1,i) + w(k,j,i) ) &
2789	* vst * rmask(j,i,sr) &
2790	* rflags_invers(j,i,k+1)
2791	ENDDO
2792	ENDDO
2793	sums_l(k,13,tn) = s1
2794	sums_l(k,15,tn) = s1
2795	ENDDO
2796	!$acc end parallel loop
2797
2798	ENDIF
2799
2800	IF ( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN
2801
2802	!$OMP DO
2803	!$acc parallel loop gang present( hom, pt, rflags_invers, rmask, sums_l, w ) create( s1 )
2804	DO k = nzb, nzt_diff
2805	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2806	DO i = nxl, nxr
2807	DO j = nys, nyn
2808	!
2809	!-- Vertical heat flux
2810	s1 = s1 + 0.5_wp * ( pt(k,j,i) - hom(k,1,4,sr) + &
2811	pt(k+1,j,i) - hom(k+1,1,4,sr) ) &
2812	* w(k,j,i) * rmask(j,i,sr) &
2813	* rflags_invers(j,i,k+1)
2814	ENDDO
2815	ENDDO
2816	sums_l(k,17,tn) = s1
2817	ENDDO
2818	!$acc end parallel loop
2819
2820	IF ( humidity ) THEN
2821
2822	!$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l, w ) create( s1 )
2823	DO k = nzb, nzt_diff
2824	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2825	DO i = nxl, nxr
2826	DO j = nys, nyn
2827	s1 = s1 + 0.5_wp * ( q(k,j,i) - hom(k,1,41,sr) + &
2828	q(k+1,j,i) - hom(k+1,1,41,sr) ) &
2829	* w(k,j,i) * rmask(j,i,sr) &
2830	* rflags_invers(j,i,k+1)
2831	ENDDO
2832	ENDDO
2833	sums_l(k,49,tn) = s1
2834	ENDDO
2835	!$acc end parallel loop
2836
2837	ENDIF
2838
2839	ENDIF
2840
2841
2842	!
2843	!-- Density at top follows Neumann condition
2844	IF ( ocean ) THEN
2845	!$acc parallel present( sums_l )
2846	sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
2847	sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
2848	!$acc end parallel
2849	ENDIF
2850
2851	!
2852	!-- Divergence of vertical flux of resolved scale energy and pressure
2853	!-- fluctuations as well as flux of pressure fluctuation itself (68).
2854	!-- First calculate the products, then the divergence.
2855	!-- Calculation is time consuming. Do it only, if profiles shall be plotted.
2856	IF ( hom(nzb+1,2,55,0) /= 0.0_wp .OR. hom(nzb+1,2,68,0) /= 0.0_wp ) THEN
2857
2858	STOP '+++ openACC porting for vertical flux div of resolved scale TKE in flow_statistics is still missing'
2859	sums_ll = 0.0_wp ! local array
2860
2861	!$OMP DO
2862	DO i = nxl, nxr
2863	DO j = nys, nyn
2864	DO k = nzb_s_inner(j,i)+1, nzt
2865
2866	sums_ll(k,1) = sums_ll(k,1) + 0.5_wp * w(k,j,i) * ( &
2867	( 0.25_wp * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1) ) &
2868	- 0.5_wp * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) ) )**2& &
2869	+ ( 0.25_wp * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i) ) &
2870	- 0.5_wp * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) ) )**2& &
2871	+ w(k,j,i)**2 )
2872
2873	sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i) &
2874	* ( p(k,j,i) + p(k+1,j,i) )
2875
2876	ENDDO
2877	ENDDO
2878	ENDDO
2879	sums_ll(0,1) = 0.0_wp ! because w is zero at the bottom
2880	sums_ll(nzt+1,1) = 0.0_wp
2881	sums_ll(0,2) = 0.0_wp
2882	sums_ll(nzt+1,2) = 0.0_wp
2883
2884	DO k = nzb+1, nzt
2885	sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
2886	sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
2887	sums_l(k,68,tn) = sums_ll(k,2)
2888	ENDDO
2889	sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
2890	sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
2891	sums_l(nzb,68,tn) = 0.0_wp ! because w* = 0 at nzb
2892
2893	ENDIF
2894
2895	!
2896	!-- Divergence of vertical flux of SGS TKE and the flux itself (69)
2897	IF ( hom(nzb+1,2,57,0) /= 0.0_wp .OR. hom(nzb+1,2,69,0) /= 0.0_wp ) THEN
2898
2899	STOP '+++ openACC porting for vertical flux div of SGS TKE in flow_statistics is still missing'
2900	!$OMP DO
2901	DO i = nxl, nxr
2902	DO j = nys, nyn
2903	DO k = nzb_s_inner(j,i)+1, nzt
2904
2905	sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5_wp * ( &
2906	(km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
2907	- (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k) &
2908	) * ddzw(k)
2909
2910	sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5_wp * ( &
2911	(km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
2912	)
2913
2914	ENDDO
2915	ENDDO
2916	ENDDO
2917	sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
2918	sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
2919
2920	ENDIF
2921
2922	!
2923	!-- Horizontal heat fluxes (subgrid, resolved, total).
2924	!-- Do it only, if profiles shall be plotted.
2925	IF ( hom(nzb+1,2,58,0) /= 0.0_wp ) THEN
2926
2927	STOP '+++ openACC porting for horizontal flux calculation in flow_statistics is still missing'
2928	!$OMP DO
2929	DO i = nxl, nxr
2930	DO j = nys, nyn
2931	DO k = nzb_s_inner(j,i)+1, nzt
2932	!
2933	!-- Subgrid horizontal heat fluxes u"pt", v"pt"
2934	sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5_wp * &
2935	( kh(k,j,i) + kh(k,j,i-1) ) &
2936	* ( pt(k,j,i-1) - pt(k,j,i) ) &
2937	* ddx * rmask(j,i,sr)
2938	sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp * &
2939	( kh(k,j,i) + kh(k,j-1,i) ) &
2940	* ( pt(k,j-1,i) - pt(k,j,i) ) &
2941	* ddy * rmask(j,i,sr)
2942	!
2943	!-- Resolved horizontal heat fluxes upt, vpt
2944	sums_l(k,59,tn) = sums_l(k,59,tn) + &
2945	( u(k,j,i) - hom(k,1,1,sr) ) * 0.5_wp * &
2946	( pt(k,j,i-1) - hom(k,1,4,sr) + &
2947	pt(k,j,i) - hom(k,1,4,sr) )
2948	pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
2949	pt(k,j,i) - hom(k,1,4,sr) )
2950	sums_l(k,62,tn) = sums_l(k,62,tn) + &
2951	( v(k,j,i) - hom(k,1,2,sr) ) * 0.5_wp * &
2952	( pt(k,j-1,i) - hom(k,1,4,sr) + &
2953	pt(k,j,i) - hom(k,1,4,sr) )
2954	ENDDO
2955	ENDDO
2956	ENDDO
2957	!
2958	!-- Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
2959	sums_l(nzb,58,tn) = 0.0_wp
2960	sums_l(nzb,59,tn) = 0.0_wp
2961	sums_l(nzb,60,tn) = 0.0_wp
2962	sums_l(nzb,61,tn) = 0.0_wp
2963	sums_l(nzb,62,tn) = 0.0_wp
2964	sums_l(nzb,63,tn) = 0.0_wp
2965
2966	ENDIF
2967
2968	!
2969	!-- Collect current large scale advection and subsidence tendencies for
2970	!-- data output
2971	IF ( large_scale_forcing .AND. ( simulated_time .GT. 0.0_wp ) ) THEN
2972	!
2973	!-- Interpolation in time of LSF_DATA
2974	nt = 1
2975	DO WHILE ( simulated_time - dt_3d > time_vert(nt) )
2976	nt = nt + 1
2977	ENDDO
2978	IF ( simulated_time - dt_3d /= time_vert(nt) ) THEN
2979	nt = nt - 1
2980	ENDIF
2981
2982	fac = ( simulated_time - dt_3d - time_vert(nt) ) &
2983	/ ( time_vert(nt+1)-time_vert(nt) )
2984
2985
2986	DO k = nzb, nzt
2987	sums_ls_l(k,0) = td_lsa_lpt(k,nt) &
2988	+ fac * ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )
2989	sums_ls_l(k,1) = td_lsa_q(k,nt) &
2990	+ fac * ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) )
2991	ENDDO
2992
2993	sums_ls_l(nzt+1,0) = sums_ls_l(nzt,0)
2994	sums_ls_l(nzt+1,1) = sums_ls_l(nzt,1)
2995
2996	IF ( large_scale_subsidence .AND. use_subsidence_tendencies ) THEN
2997
2998	DO k = nzb, nzt
2999	sums_ls_l(k,2) = td_sub_lpt(k,nt) + fac * &
3000	( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )
3001	sums_ls_l(k,3) = td_sub_q(k,nt) + fac * &
3002	( td_sub_q(k,nt+1) - td_sub_q(k,nt) )
3003	ENDDO
3004
3005	sums_ls_l(nzt+1,2) = sums_ls_l(nzt,2)
3006	sums_ls_l(nzt+1,3) = sums_ls_l(nzt,3)
3007
3008	ENDIF
3009
3010	ENDIF
3011
3012
3013	IF ( land_surface ) THEN
3014	!$OMP DO
3015	DO i = nxl, nxr
3016	DO j = nys, nyn
3017	DO k = nzb_soil, nzt_soil
3018	sums_l(k,89,tn) = sums_l(k,89,tn) + t_soil(k,j,i) * rmask(j,i,sr)
3019	sums_l(k,91,tn) = sums_l(k,91,tn) + m_soil(k,j,i) * rmask(j,i,sr)
3020	ENDDO
3021	ENDDO
3022	ENDDO
3023	ENDIF
3024
3025
3026	IF ( radiation .AND. radiation_scheme == 'rrtmg' ) THEN
3027	!$OMP DO
3028	DO i = nxl, nxr
3029	DO j = nys, nyn
3030	DO k = nzb_s_inner(j,i)+1, nzt+1
3031	sums_l(k,102,tn) = sums_l(k,102,tn) + rad_lw_in(k,j,i) * rmask(j,i,sr)
3032	sums_l(k,103,tn) = sums_l(k,103,tn) + rad_lw_out(k,j,i) * rmask(j,i,sr)
3033	sums_l(k,104,tn) = sums_l(k,104,tn) + rad_sw_in(k,j,i) * rmask(j,i,sr)
3034	sums_l(k,105,tn) = sums_l(k,105,tn) + rad_sw_out(k,j,i) * rmask(j,i,sr)
3035	ENDDO
3036	ENDDO
3037	ENDDO
3038	ENDIF
3039
3040	!
3041	!-- Calculate the user-defined profiles
3042	CALL user_statistics( 'profiles', sr, tn )
3043	!$OMP END PARALLEL
3044
3045	!
3046	!-- Summation of thread sums
3047	IF ( threads_per_task > 1 ) THEN
3048	STOP '+++ openACC porting for threads_per_task > 1 in flow_statistics is still missing'
3049	DO i = 1, threads_per_task-1
3050	sums_l(:,3,0) = sums_l(:,3,0) + sums_l(:,3,i)
3051	sums_l(:,4:40,0) = sums_l(:,4:40,0) + sums_l(:,4:40,i)
3052	sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
3053	sums_l(:,45:pr_palm,i)
3054	IF ( max_pr_user > 0 ) THEN
3055	sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
3056	sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
3057	sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
3058	ENDIF
3059	ENDDO
3060	ENDIF
3061
3062	!$acc update host( hom, sums, sums_l )
3063
3064	#if defined( __parallel )
3065
3066	!
3067	!-- Compute total sum from local sums
3068	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
3069	CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, &
3070	MPI_SUM, comm2d, ierr )
3071	IF ( large_scale_forcing ) THEN
3072	CALL MPI_ALLREDUCE( sums_ls_l(nzb,2), sums(nzb,83), ngp_sums_ls, &
3073	MPI_REAL, MPI_SUM, comm2d, ierr )
3074	ENDIF
3075	#else
3076	sums = sums_l(:,:,0)
3077	IF ( large_scale_forcing ) THEN
3078	sums(:,81:88) = sums_ls_l
3079	ENDIF
3080	#endif
3081
3082	!
3083	!-- Final values are obtained by division by the total number of grid points
3084	!-- used for summation. After that store profiles.
3085	!-- Profiles:
3086	DO k = nzb, nzt+1
3087	sums(k,3) = sums(k,3) / ngp_2dh(sr)
3088	sums(k,8:11) = sums(k,8:11) / ngp_2dh_s_inner(k,sr)
3089	sums(k,12:22) = sums(k,12:22) / ngp_2dh(sr)
3090	sums(k,23:29) = sums(k,23:29) / ngp_2dh_s_inner(k,sr)
3091	sums(k,30:32) = sums(k,30:32) / ngp_2dh(sr)
3092	sums(k,33:34) = sums(k,33:34) / ngp_2dh_s_inner(k,sr)
3093	sums(k,35:39) = sums(k,35:39) / ngp_2dh(sr)
3094	sums(k,40) = sums(k,40) / ngp_2dh_s_inner(k,sr)
3095	sums(k,45:53) = sums(k,45:53) / ngp_2dh(sr)
3096	sums(k,54) = sums(k,54) / ngp_2dh_s_inner(k,sr)
3097	sums(k,55:63) = sums(k,55:63) / ngp_2dh(sr)
3098	sums(k,64) = sums(k,64) / ngp_2dh_s_inner(k,sr)
3099	sums(k,70:80) = sums(k,70:80) / ngp_2dh_s_inner(k,sr)
3100	sums(k,81:88) = sums(k,81:88) / ngp_2dh(sr)
3101	sums(k,89:105) = sums(k,89:105) / ngp_2dh(sr)
3102	sums(k,106:pr_palm-2) = sums(k,106:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
3103	ENDDO
3104
3105	!-- Upstream-parts
3106	sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
3107	!-- u* and so on
3108	!-- As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
3109	!-- size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
3110	!-- above the topography, they are being divided by ngp_2dh(sr)
3111	sums(nzb:nzb+3,pr_palm) = sums(nzb:nzb+3,pr_palm) / &
3112	ngp_2dh(sr)
3113	sums(nzb+12,pr_palm) = sums(nzb+12,pr_palm) / & ! qs
3114	ngp_2dh(sr)
3115	!-- eges, e*
3116	sums(nzb+4:nzb+5,pr_palm) = sums(nzb+4:nzb+5,pr_palm) / &
3117	ngp_3d(sr)
3118	!-- Old and new divergence
3119	sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
3120	ngp_3d_inner(sr)
3121
3122	!-- User-defined profiles
3123	IF ( max_pr_user > 0 ) THEN
3124	DO k = nzb, nzt+1
3125	sums(k,pr_palm+1:pr_palm+max_pr_user) = &
3126	sums(k,pr_palm+1:pr_palm+max_pr_user) / &
3127	ngp_2dh_s_inner(k,sr)
3128	ENDDO
3129	ENDIF
3130
3131	!
3132	!-- Collect horizontal average in hom.
3133	!-- Compute deduced averages (e.g. total heat flux)
3134	hom(:,1,3,sr) = sums(:,3) ! w
3135	hom(:,1,8,sr) = sums(:,8) ! e profiles 5-7 are initial profiles
3136	hom(:,1,9,sr) = sums(:,9) ! km
3137	hom(:,1,10,sr) = sums(:,10) ! kh
3138	hom(:,1,11,sr) = sums(:,11) ! l
3139	hom(:,1,12,sr) = sums(:,12) ! w"u"
3140	hom(:,1,13,sr) = sums(:,13) ! wu
3141	hom(:,1,14,sr) = sums(:,14) ! w"v"
3142	hom(:,1,15,sr) = sums(:,15) ! wv
3143	hom(:,1,16,sr) = sums(:,16) ! w"pt"
3144	hom(:,1,17,sr) = sums(:,17) ! wpt
3145	hom(:,1,18,sr) = sums(:,16) + sums(:,17) ! wpt
3146	hom(:,1,19,sr) = sums(:,12) + sums(:,13) ! wu
3147	hom(:,1,20,sr) = sums(:,14) + sums(:,15) ! wv
3148	hom(:,1,21,sr) = sums(:,21) ! wptBC
3149	hom(:,1,22,sr) = sums(:,16) + sums(:,21) ! wptBC
3150	! profile 24 is initial profile (sa)
3151	! profiles 25-29 left empty for initial
3152	! profiles
3153	hom(:,1,30,sr) = sums(:,30) ! u*2
3154	hom(:,1,31,sr) = sums(:,31) ! v*2
3155	hom(:,1,32,sr) = sums(:,32) ! w*2
3156	hom(:,1,33,sr) = sums(:,33) ! pt*2
3157	hom(:,1,34,sr) = sums(:,34) ! e*
3158	hom(:,1,35,sr) = sums(:,35) ! w2pt
3159	hom(:,1,36,sr) = sums(:,36) ! wpt2
3160	hom(:,1,37,sr) = sums(:,37) ! we
3161	hom(:,1,38,sr) = sums(:,38) ! w*3
3162	hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20_wp )**1.5_wp ! Sw
3163	hom(:,1,40,sr) = sums(:,40) ! p
3164	hom(:,1,45,sr) = sums(:,45) ! w"vpt"
3165	hom(:,1,46,sr) = sums(:,46) ! wvpt
3166	hom(:,1,47,sr) = sums(:,45) + sums(:,46) ! wvpt
3167	hom(:,1,48,sr) = sums(:,48) ! w"q" (w"qv")
3168	hom(:,1,49,sr) = sums(:,49) ! wq (wqv)
3169	hom(:,1,50,sr) = sums(:,48) + sums(:,49) ! wq (wqv)
3170	hom(:,1,51,sr) = sums(:,51) ! w"qv"
3171	hom(:,1,52,sr) = sums(:,52) ! wqv
3172	hom(:,1,53,sr) = sums(:,52) + sums(:,51) ! wq (wqv)
3173	hom(:,1,54,sr) = sums(:,54) ! ql
3174	hom(:,1,55,sr) = sums(:,55) ! wuu*/dz
3175	hom(:,1,56,sr) = sums(:,56) ! wp/dz
3176	hom(:,1,57,sr) = sums(:,57) ! ( w"e + w"p"/rho )/dz
3177	hom(:,1,58,sr) = sums(:,58) ! u"pt"
3178	hom(:,1,59,sr) = sums(:,59) ! upt
3179	hom(:,1,60,sr) = sums(:,58) + sums(:,59) ! upt_t
3180	hom(:,1,61,sr) = sums(:,61) ! v"pt"
3181	hom(:,1,62,sr) = sums(:,62) ! vpt
3182	hom(:,1,63,sr) = sums(:,61) + sums(:,62) ! vpt_t
3183	hom(:,1,64,sr) = sums(:,64) ! rho
3184	hom(:,1,65,sr) = sums(:,65) ! w"sa"
3185	hom(:,1,66,sr) = sums(:,66) ! wsa
3186	hom(:,1,67,sr) = sums(:,65) + sums(:,66) ! wsa
3187	hom(:,1,68,sr) = sums(:,68) ! wp
3188	hom(:,1,69,sr) = sums(:,69) ! w"e + w"p"/rho
3189	hom(:,1,70,sr) = sums(:,70) ! q*2
3190	hom(:,1,71,sr) = sums(:,71) ! prho
3191	hom(:,1,72,sr) = hyp * 1E-4_wp ! hyp in dbar
3192	hom(:,1,73,sr) = sums(:,73) ! nr
3193	hom(:,1,74,sr) = sums(:,74) ! qr
3194	hom(:,1,75,sr) = sums(:,75) ! qc
3195	hom(:,1,76,sr) = sums(:,76) ! prr (precipitation rate)
3196	! 77 is initial density profile
3197	hom(:,1,78,sr) = ug ! ug
3198	hom(:,1,79,sr) = vg ! vg
3199	hom(:,1,80,sr) = w_subs ! w_subs
3200
3201	IF ( large_scale_forcing ) THEN
3202	hom(:,1,81,sr) = sums_ls_l(:,0) ! td_lsa_lpt
3203	hom(:,1,82,sr) = sums_ls_l(:,1) ! td_lsa_q
3204	IF ( use_subsidence_tendencies ) THEN
3205	hom(:,1,83,sr) = sums_ls_l(:,2) ! td_sub_lpt
3206	hom(:,1,84,sr) = sums_ls_l(:,3) ! td_sub_q
3207	ELSE
3208	hom(:,1,83,sr) = sums(:,83) ! td_sub_lpt
3209	hom(:,1,84,sr) = sums(:,84) ! td_sub_q
3210	ENDIF
3211	hom(:,1,85,sr) = sums(:,85) ! td_nud_lpt
3212	hom(:,1,86,sr) = sums(:,86) ! td_nud_q
3213	hom(:,1,87,sr) = sums(:,87) ! td_nud_u
3214	hom(:,1,88,sr) = sums(:,88) ! td_nud_v
3215	END IF
3216
3217	hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
3218	! upstream-parts u_x, u_y, u_z, v_x,
3219	! v_y, usw. (in last but one profile)
3220	hom(:,1,pr_palm,sr) = sums(:,pr_palm)
3221	! u, w'u', w'v', t (in last profile)
3222
3223	IF ( max_pr_user > 0 ) THEN ! user-defined profiles
3224	hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
3225	sums(:,pr_palm+1:pr_palm+max_pr_user)
3226	ENDIF
3227
3228	!
3229	!-- Determine the boundary layer height using two different schemes.
3230	!-- First scheme: Starting from the Earth's (Ocean's) surface, look for the
3231	!-- first relative minimum (maximum) of the total heat flux.
3232	!-- The corresponding height is assumed as the boundary layer height, if it
3233	!-- is less than 1.5 times the height where the heat flux becomes negative
3234	!-- (positive) for the first time.
3235	z_i(1) = 0.0_wp
3236	first = .TRUE.
3237
3238	IF ( ocean ) THEN
3239	DO k = nzt, nzb+1, -1
3240	IF ( first .AND. hom(k,1,18,sr) < 0.0_wp &
3241	.AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp ) THEN
3242	first = .FALSE.
3243	height = zw(k)
3244	ENDIF
3245	IF ( hom(k,1,18,sr) < 0.0_wp .AND. &
3246	abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND. &
3247	hom(k-1,1,18,sr) > hom(k,1,18,sr) ) THEN
3248	IF ( zw(k) < 1.5_wp * height ) THEN
3249	z_i(1) = zw(k)
3250	ELSE
3251	z_i(1) = height
3252	ENDIF
3253	EXIT
3254	ENDIF
3255	ENDDO
3256	ELSE
3257	DO k = nzb, nzt-1
3258	IF ( first .AND. hom(k,1,18,sr) < 0.0_wp &
3259	.AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp ) THEN
3260	first = .FALSE.
3261	height = zw(k)
3262	ENDIF
3263	IF ( hom(k,1,18,sr) < 0.0 .AND. &
3264	abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND. &
3265	hom(k+1,1,18,sr) > hom(k,1,18,sr) ) THEN
3266	IF ( zw(k) < 1.5_wp * height ) THEN
3267	z_i(1) = zw(k)
3268	ELSE
3269	z_i(1) = height
3270	ENDIF
3271	EXIT
3272	ENDIF
3273	ENDDO
3274	ENDIF
3275
3276	!
3277	!-- Second scheme: Gradient scheme from Sullivan et al. (1998), modified
3278	!-- by Uhlenbrock(2006). The boundary layer height is the height with the
3279	!-- maximal local temperature gradient: starting from the second (the last
3280	!-- but one) vertical gridpoint, the local gradient must be at least
3281	!-- 0.2K/100m and greater than the next four gradients.
3282	!-- WARNING: The threshold value of 0.2K/100m must be adjusted for the
3283	!-- ocean case!
3284	z_i(2) = 0.0_wp
3285	DO k = nzb+1, nzt+1
3286	dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
3287	ENDDO
3288	dptdz_threshold = 0.2_wp / 100.0_wp
3289
3290	IF ( ocean ) THEN
3291	DO k = nzt+1, nzb+5, -1
3292	IF ( dptdz(k) > dptdz_threshold .AND. &
3293	dptdz(k) > dptdz(k-1) .AND. dptdz(k) > dptdz(k-2) .AND. &
3294	dptdz(k) > dptdz(k-3) .AND. dptdz(k) > dptdz(k-4) ) THEN
3295	z_i(2) = zw(k-1)
3296	EXIT
3297	ENDIF
3298	ENDDO
3299	ELSE
3300	DO k = nzb+1, nzt-3
3301	IF ( dptdz(k) > dptdz_threshold .AND. &
3302	dptdz(k) > dptdz(k+1) .AND. dptdz(k) > dptdz(k+2) .AND. &
3303	dptdz(k) > dptdz(k+3) .AND. dptdz(k) > dptdz(k+4) ) THEN
3304	z_i(2) = zw(k-1)
3305	EXIT
3306	ENDIF
3307	ENDDO
3308	ENDIF
3309
3310	hom(nzb+6,1,pr_palm,sr) = z_i(1)
3311	hom(nzb+7,1,pr_palm,sr) = z_i(2)
3312
3313	!
3314	!-- Computation of both the characteristic vertical velocity and
3315	!-- the characteristic convective boundary layer temperature.
3316	!-- The horizontal average at nzb+1 is input for the average temperature.
3317	IF ( hom(nzb,1,18,sr) > 0.0_wp .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8_wp &
3318	.AND. z_i(1) /= 0.0_wp ) THEN
3319	hom(nzb+8,1,pr_palm,sr) = ( g / hom(nzb+1,1,4,sr) * &
3320	hom(nzb,1,18,sr) * &
3321	ABS( z_i(1) ) )**0.333333333_wp
3322	!-- so far this only works if Prandtl layer is used
3323	hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
3324	ELSE
3325	hom(nzb+8,1,pr_palm,sr) = 0.0_wp
3326	hom(nzb+11,1,pr_palm,sr) = 0.0_wp
3327	ENDIF
3328
3329	!
3330	!-- Collect the time series quantities
3331	ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr) ! E
3332	ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr) ! E*
3333	ts_value(3,sr) = dt_3d
3334	ts_value(4,sr) = hom(nzb,1,pr_palm,sr) ! u*
3335	ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr) ! th*
3336	ts_value(6,sr) = u_max
3337	ts_value(7,sr) = v_max
3338	ts_value(8,sr) = w_max
3339	ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr) ! new divergence
3340	ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr) ! old Divergence
3341	ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr) ! z_i(1)
3342	ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr) ! z_i(2)
3343	ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr) ! w*
3344	ts_value(14,sr) = hom(nzb,1,16,sr) ! w'pt' at k=0
3345	ts_value(15,sr) = hom(nzb+1,1,16,sr) ! w'pt' at k=1
3346	ts_value(16,sr) = hom(nzb+1,1,18,sr) ! wpt at k=1
3347	ts_value(17,sr) = hom(nzb,1,4,sr) ! pt(0)
3348	ts_value(18,sr) = hom(nzb+1,1,4,sr) ! pt(zp)
3349	ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr) ! u'w' at k=0
3350	ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr) ! v'w' at k=0
3351	ts_value(21,sr) = hom(nzb,1,48,sr) ! w"q" at k=0
3352
3353	IF ( ts_value(5,sr) /= 0.0_wp ) THEN
3354	ts_value(22,sr) = ts_value(4,sr)**2_wp / &
3355	( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
3356	ELSE
3357	ts_value(22,sr) = 10000.0_wp
3358	ENDIF
3359
3360	ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr) ! q*
3361
3362	!
3363	!-- Collect land surface model timeseries
3364	IF ( land_surface ) THEN
3365	ts_value(dots_soil ,sr) = hom(nzb,1,93,sr) ! ghf_eb
3366	ts_value(dots_soil+1,sr) = hom(nzb,1,94,sr) ! shf_eb
3367	ts_value(dots_soil+2,sr) = hom(nzb,1,95,sr) ! qsws_eb
3368	ts_value(dots_soil+3,sr) = hom(nzb,1,96,sr) ! qsws_liq_eb
3369	ts_value(dots_soil+4,sr) = hom(nzb,1,97,sr) ! qsws_soil_eb
3370	ts_value(dots_soil+5,sr) = hom(nzb,1,98,sr) ! qsws_veg_eb
3371	ts_value(dots_soil+6,sr) = hom(nzb,1,99,sr) ! r_a
3372	ts_value(dots_soil+7,sr) = hom(nzb,1,100,sr) ! r_s
3373	ENDIF
3374	!
3375	!-- Collect radiation model timeseries
3376	IF ( radiation ) THEN
3377	ts_value(dots_rad,sr) = hom(nzb,1,101,sr) ! rad_net
3378	ts_value(dots_rad+1,sr) = hom(nzb,1,102,sr) ! rad_lw_in
3379	ts_value(dots_rad+2,sr) = hom(nzb,1,103,sr) ! rad_lw_out
3380	ts_value(dots_rad+3,sr) = hom(nzb,1,104,sr) ! rad_lw_in
3381	ts_value(dots_rad+4,sr) = hom(nzb,1,105,sr) ! rad_lw_out
3382
3383	#if defined ( __rrtmg )
3384	IF ( radiation_scheme == 'rrtmg' ) THEN
3385	ts_value(dots_rad+5,sr) = hom(nzb,1,106,sr) ! rrtm_aldif
3386	ts_value(dots_rad+6,sr) = hom(nzb,1,107,sr) ! rrtm_aldir
3387	ts_value(dots_rad+7,sr) = hom(nzb,1,108,sr) ! rrtm_asdif
3388	ts_value(dots_rad+8,sr) = hom(nzb,1,109,sr) ! rrtm_asdir
3389	ENDIF
3390	#endif
3391
3392	ENDIF
3393
3394	!
3395	!-- Calculate additional statistics provided by the user interface
3396	CALL user_statistics( 'time_series', sr, 0 )
3397
3398	ENDDO ! loop of the subregions
3399
3400	!$acc end data
3401
3402	!
3403	!-- If required, sum up horizontal averages for subsequent time averaging
3404	IF ( do_sum ) THEN
3405	IF ( average_count_pr == 0 ) hom_sum = 0.0_wp
3406	hom_sum = hom_sum + hom(:,1,:,:)
3407	average_count_pr = average_count_pr + 1
3408	do_sum = .FALSE.
3409	ENDIF
3410
3411	!
3412	!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
3413	!-- This flag is reset after each time step in time_integration.
3414	flow_statistics_called = .TRUE.
3415
3416	CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
3417
3418
3419	END SUBROUTINE flow_statistics
3420	#endif

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |