Home

Context Navigation

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

Last change on this file since 1463 was 1451, checked in by heinze, 10 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 131.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |