Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |