Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |