Home

Context Navigation

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

Last change on this file since 1353 was 1353, checked in by heinze, 10 years ago
REAL constants provided with KIND-attribute
Property svn:keywords set to `Id`
File size: 124.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |