Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |