Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |