Home

Context Navigation

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

Last change on this file since 1386 was 1386, checked in by boeske, 10 years ago
bugfix: simulated time before last timestep must be used for the profile data output of large scale forcing tendencies
Property svn:keywords set to `Id`
File size: 130.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |