Home

Context Navigation

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

Last change on this file since 1330 was 1323, checked in by raasch, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 121.3 KB

Line
1	#if ! defined( __openacc )
2	SUBROUTINE flow_statistics
3
4	!--------------------------------------------------------------------------------!
5	! This file is part of PALM.
6	!
7	! PALM is free software: you can redistribute it and/or modify it under the terms
8	! of the GNU General Public License as published by the Free Software Foundation,
9	! either version 3 of the License, or (at your option) any later version.
10	!
11	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
12	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
13	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
14	!
15	! You should have received a copy of the GNU General Public License along with
16	! PALM. If not, see <http://www.gnu.org/licenses/>.
17	!
18	! Copyright 1997-2014 Leibniz Universitaet Hannover
19	!--------------------------------------------------------------------------------!
20	!
21	! Current revisions:
22	! -----------------
23	!
24	!
25	! Former revisions:
26	! -----------------
27	! $Id: flow_statistics.f90 1323 2014-03-20 17:09:54Z suehring $
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
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 * &
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 * ( &
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
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 * ( 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
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 * &
520	( ust2 + vst2 + w2 ) * rmask(j,i,sr)
521	sums_l_eper = sums_l_eper + &
522	0.5 * ( 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 * ( &
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 * ( &
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 * ( 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 * ( 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 * ( 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 * ( 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 * ( 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 * ( 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 * rmask(j,i,sr) ! u"pt"
623	sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
624	0.0 * 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 + 0.61 * q(nzb,j,i) ) * &
634	shf(j,i) + 0.61 * 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 + 0.61 * q(nzb,j,i) - &
639	ql(nzb,j,i) ) * shf(j,i) + &
640	0.61 * 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 * rmask(j,i,sr) ! u"pt"
666	sums_l(nzt:nzt+1,61,tn) = sums_l(nzt:nzt+1,61,tn) + &
667	0.0 * 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 + 0.61 * q(nzt,j,i) ) * &
678	tswst(j,i) + 0.61 * 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 + 0.61 * q(nzt,j,i) - &
683	ql(nzt,j,i) ) * tswst(j,i) + &
684	0.61 * 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 * ( 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 * ( 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 * ( 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 * ( 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 * ( 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 * &
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 * ( 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 + 0.61 * hom(k,1,41,sr) ) * &
780	sums_l(k,17,tn) + &
781	0.61 * hom(k,1,4,sr) * sums_l(k,49,tn)
782	END IF
783	END IF
784	ENDIF
785	!
786	!-- Passive scalar flux
787	IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca &
788	.OR. sr /= 0 ) ) THEN
789	pts = 0.5 * ( q(k,j,i) - hom(k,1,41,sr) + &
790	q(k+1,j,i) - hom(k+1,1,41,sr) )
791	sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
792	rmask(j,i,sr)
793	ENDIF
794
795	!
796	!-- Energy flux we
797	!-- has to be adjusted
798	sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5 * &
799	( ust2 + vst2 + w(k,j,i)**2 )&
800	* rmask(j,i,sr)
801	ENDDO
802	ENDDO
803	ENDDO
804	!
805	!-- For speed optimization fluxes which have been computed in part directly
806	!-- inside the WS advection routines are treated seperatly
807	!-- Momentum fluxes first:
808	IF ( .NOT. ws_scheme_mom .OR. sr /= 0 ) THEN
809	!$OMP DO
810	DO i = nxl, nxr
811	DO j = nys, nyn
812	DO k = nzb_diff_s_inner(j,i)-1, nzt_diff
813	ust = 0.5 * ( u(k,j,i) - hom(k,1,1,sr) + &
814	u(k+1,j,i) - hom(k+1,1,1,sr) )
815	vst = 0.5 * ( v(k,j,i) - hom(k,1,2,sr) + &
816	v(k+1,j,i) - hom(k+1,1,2,sr) )
817	!
818	!-- Momentum flux wu
819	sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5 * &
820	( w(k,j,i-1) + w(k,j,i) ) &
821	* ust * rmask(j,i,sr)
822	!
823	!-- Momentum flux wv
824	sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5 * &
825	( w(k,j-1,i) + w(k,j,i) ) &
826	* vst * rmask(j,i,sr)
827	ENDDO
828	ENDDO
829	ENDDO
830
831	ENDIF
832	IF ( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN
833	!$OMP DO
834	DO i = nxl, nxr
835	DO j = nys, nyn
836	DO k = nzb_diff_s_inner(j,i)-1, nzt_diff
837	!
838	!-- Vertical heat flux
839	sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5 * &
840	( pt(k,j,i) - hom(k,1,4,sr) + &
841	pt(k+1,j,i) - hom(k+1,1,4,sr) ) &
842	* w(k,j,i) * rmask(j,i,sr)
843	IF ( humidity ) THEN
844	pts = 0.5 * ( q(k,j,i) - hom(k,1,41,sr) + &
845	q(k+1,j,i) - hom(k+1,1,41,sr) )
846	sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
847	rmask(j,i,sr)
848	ENDIF
849	ENDDO
850	ENDDO
851	ENDDO
852
853	ENDIF
854
855	!
856	!-- Density at top follows Neumann condition
857	IF ( ocean ) THEN
858	sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
859	sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
860	ENDIF
861
862	!
863	!-- Divergence of vertical flux of resolved scale energy and pressure
864	!-- fluctuations as well as flux of pressure fluctuation itself (68).
865	!-- First calculate the products, then the divergence.
866	!-- Calculation is time consuming. Do it only, if profiles shall be plotted.
867	IF ( hom(nzb+1,2,55,0) /= 0.0 .OR. hom(nzb+1,2,68,0) /= 0.0 ) THEN
868
869	sums_ll = 0.0 ! local array
870
871	!$OMP DO
872	DO i = nxl, nxr
873	DO j = nys, nyn
874	DO k = nzb_s_inner(j,i)+1, nzt
875
876	sums_ll(k,1) = sums_ll(k,1) + 0.5 * w(k,j,i) * ( &
877	( 0.25 * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1) &
878	- 0.5 * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) ) &
879	) )**2 &
880	+ ( 0.25 * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i) &
881	- 0.5 * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) ) &
882	) )**2 &
883	+ w(k,j,i)**2 )
884
885	sums_ll(k,2) = sums_ll(k,2) + 0.5 * w(k,j,i) &
886	* ( p(k,j,i) + p(k+1,j,i) )
887
888	ENDDO
889	ENDDO
890	ENDDO
891	sums_ll(0,1) = 0.0 ! because w is zero at the bottom
892	sums_ll(nzt+1,1) = 0.0
893	sums_ll(0,2) = 0.0
894	sums_ll(nzt+1,2) = 0.0
895
896	DO k = nzb+1, nzt
897	sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
898	sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
899	sums_l(k,68,tn) = sums_ll(k,2)
900	ENDDO
901	sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
902	sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
903	sums_l(nzb,68,tn) = 0.0 ! because w* = 0 at nzb
904
905	ENDIF
906
907	!
908	!-- Divergence of vertical flux of SGS TKE and the flux itself (69)
909	IF ( hom(nzb+1,2,57,0) /= 0.0 .OR. hom(nzb+1,2,69,0) /= 0.0 ) THEN
910
911	!$OMP DO
912	DO i = nxl, nxr
913	DO j = nys, nyn
914	DO k = nzb_s_inner(j,i)+1, nzt
915
916	sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5 * ( &
917	(km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
918	- (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k) &
919	) * ddzw(k)
920
921	sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5 * ( &
922	(km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
923	)
924
925	ENDDO
926	ENDDO
927	ENDDO
928	sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
929	sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
930
931	ENDIF
932
933	!
934	!-- Horizontal heat fluxes (subgrid, resolved, total).
935	!-- Do it only, if profiles shall be plotted.
936	IF ( hom(nzb+1,2,58,0) /= 0.0 ) THEN
937
938	!$OMP DO
939	DO i = nxl, nxr
940	DO j = nys, nyn
941	DO k = nzb_s_inner(j,i)+1, nzt
942	!
943	!-- Subgrid horizontal heat fluxes u"pt", v"pt"
944	sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5 * &
945	( kh(k,j,i) + kh(k,j,i-1) ) &
946	* ( pt(k,j,i-1) - pt(k,j,i) ) &
947	* ddx * rmask(j,i,sr)
948	sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5 * &
949	( kh(k,j,i) + kh(k,j-1,i) ) &
950	* ( pt(k,j-1,i) - pt(k,j,i) ) &
951	* ddy * rmask(j,i,sr)
952	!
953	!-- Resolved horizontal heat fluxes upt, vpt
954	sums_l(k,59,tn) = sums_l(k,59,tn) + &
955	( u(k,j,i) - hom(k,1,1,sr) ) &
956	* 0.5 * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
957	pt(k,j,i) - hom(k,1,4,sr) )
958	pts = 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
959	pt(k,j,i) - hom(k,1,4,sr) )
960	sums_l(k,62,tn) = sums_l(k,62,tn) + &
961	( v(k,j,i) - hom(k,1,2,sr) ) &
962	* 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
963	pt(k,j,i) - hom(k,1,4,sr) )
964	ENDDO
965	ENDDO
966	ENDDO
967	!
968	!-- Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
969	sums_l(nzb,58,tn) = 0.0
970	sums_l(nzb,59,tn) = 0.0
971	sums_l(nzb,60,tn) = 0.0
972	sums_l(nzb,61,tn) = 0.0
973	sums_l(nzb,62,tn) = 0.0
974	sums_l(nzb,63,tn) = 0.0
975
976	ENDIF
977
978	!
979	!-- Calculate the user-defined profiles
980	CALL user_statistics( 'profiles', sr, tn )
981	!$OMP END PARALLEL
982
983	!
984	!-- Summation of thread sums
985	IF ( threads_per_task > 1 ) THEN
986	DO i = 1, threads_per_task-1
987	sums_l(:,3,0) = sums_l(:,3,0) + sums_l(:,3,i)
988	sums_l(:,4:40,0) = sums_l(:,4:40,0) + sums_l(:,4:40,i)
989	sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
990	sums_l(:,45:pr_palm,i)
991	IF ( max_pr_user > 0 ) THEN
992	sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
993	sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
994	sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
995	ENDIF
996	ENDDO
997	ENDIF
998
999	#if defined( __parallel )
1000
1001	!
1002	!-- Compute total sum from local sums
1003	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1004	CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, &
1005	MPI_SUM, comm2d, ierr )
1006	#else
1007	sums = sums_l(:,:,0)
1008	#endif
1009
1010	!
1011	!-- Final values are obtained by division by the total number of grid points
1012	!-- used for summation. After that store profiles.
1013	!-- Profiles:
1014	DO k = nzb, nzt+1
1015	sums(k,3) = sums(k,3) / ngp_2dh(sr)
1016	sums(k,8:11) = sums(k,8:11) / ngp_2dh_s_inner(k,sr)
1017	sums(k,12:22) = sums(k,12:22) / ngp_2dh(sr)
1018	sums(k,23:29) = sums(k,23:29) / ngp_2dh_s_inner(k,sr)
1019	sums(k,30:32) = sums(k,30:32) / ngp_2dh(sr)
1020	sums(k,33:34) = sums(k,33:34) / ngp_2dh_s_inner(k,sr)
1021	sums(k,35:39) = sums(k,35:39) / ngp_2dh(sr)
1022	sums(k,40) = sums(k,40) / ngp_2dh_s_inner(k,sr)
1023	sums(k,45:53) = sums(k,45:53) / ngp_2dh(sr)
1024	sums(k,54) = sums(k,54) / ngp_2dh_s_inner(k,sr)
1025	sums(k,55:63) = sums(k,55:63) / ngp_2dh(sr)
1026	sums(k,64) = sums(k,64) / ngp_2dh_s_inner(k,sr)
1027	sums(k,65:69) = sums(k,65:69) / ngp_2dh(sr)
1028	sums(k,70:pr_palm-2) = sums(k,70:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
1029	ENDDO
1030
1031	!-- Upstream-parts
1032	sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
1033	!-- u* and so on
1034	!-- As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
1035	!-- size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
1036	!-- above the topography, they are being divided by ngp_2dh(sr)
1037	sums(nzb:nzb+3,pr_palm) = sums(nzb:nzb+3,pr_palm) / &
1038	ngp_2dh(sr)
1039	sums(nzb+12,pr_palm) = sums(nzb+12,pr_palm) / & ! qs
1040	ngp_2dh(sr)
1041	!-- eges, e*
1042	sums(nzb+4:nzb+5,pr_palm) = sums(nzb+4:nzb+5,pr_palm) / &
1043	ngp_3d(sr)
1044	!-- Old and new divergence
1045	sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
1046	ngp_3d_inner(sr)
1047
1048	!-- User-defined profiles
1049	IF ( max_pr_user > 0 ) THEN
1050	DO k = nzb, nzt+1
1051	sums(k,pr_palm+1:pr_palm+max_pr_user) = &
1052	sums(k,pr_palm+1:pr_palm+max_pr_user) / &
1053	ngp_2dh_s_inner(k,sr)
1054	ENDDO
1055	ENDIF
1056
1057	!
1058	!-- Collect horizontal average in hom.
1059	!-- Compute deduced averages (e.g. total heat flux)
1060	hom(:,1,3,sr) = sums(:,3) ! w
1061	hom(:,1,8,sr) = sums(:,8) ! e profiles 5-7 are initial profiles
1062	hom(:,1,9,sr) = sums(:,9) ! km
1063	hom(:,1,10,sr) = sums(:,10) ! kh
1064	hom(:,1,11,sr) = sums(:,11) ! l
1065	hom(:,1,12,sr) = sums(:,12) ! w"u"
1066	hom(:,1,13,sr) = sums(:,13) ! wu
1067	hom(:,1,14,sr) = sums(:,14) ! w"v"
1068	hom(:,1,15,sr) = sums(:,15) ! wv
1069	hom(:,1,16,sr) = sums(:,16) ! w"pt"
1070	hom(:,1,17,sr) = sums(:,17) ! wpt
1071	hom(:,1,18,sr) = sums(:,16) + sums(:,17) ! wpt
1072	hom(:,1,19,sr) = sums(:,12) + sums(:,13) ! wu
1073	hom(:,1,20,sr) = sums(:,14) + sums(:,15) ! wv
1074	hom(:,1,21,sr) = sums(:,21) ! wptBC
1075	hom(:,1,22,sr) = sums(:,16) + sums(:,21) ! wptBC
1076	! profile 24 is initial profile (sa)
1077	! profiles 25-29 left empty for initial
1078	! profiles
1079	hom(:,1,30,sr) = sums(:,30) ! u*2
1080	hom(:,1,31,sr) = sums(:,31) ! v*2
1081	hom(:,1,32,sr) = sums(:,32) ! w*2
1082	hom(:,1,33,sr) = sums(:,33) ! pt*2
1083	hom(:,1,34,sr) = sums(:,34) ! e*
1084	hom(:,1,35,sr) = sums(:,35) ! w2pt
1085	hom(:,1,36,sr) = sums(:,36) ! wpt2
1086	hom(:,1,37,sr) = sums(:,37) ! we
1087	hom(:,1,38,sr) = sums(:,38) ! w*3
1088	hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20 )**1.5 ! Sw
1089	hom(:,1,40,sr) = sums(:,40) ! p
1090	hom(:,1,45,sr) = sums(:,45) ! w"vpt"
1091	hom(:,1,46,sr) = sums(:,46) ! wvpt
1092	hom(:,1,47,sr) = sums(:,45) + sums(:,46) ! wvpt
1093	hom(:,1,48,sr) = sums(:,48) ! w"q" (w"qv")
1094	hom(:,1,49,sr) = sums(:,49) ! wq (wqv)
1095	hom(:,1,50,sr) = sums(:,48) + sums(:,49) ! wq (wqv)
1096	hom(:,1,51,sr) = sums(:,51) ! w"qv"
1097	hom(:,1,52,sr) = sums(:,52) ! wqv
1098	hom(:,1,53,sr) = sums(:,52) + sums(:,51) ! wq (wqv)
1099	hom(:,1,54,sr) = sums(:,54) ! ql
1100	hom(:,1,55,sr) = sums(:,55) ! wuu*/dz
1101	hom(:,1,56,sr) = sums(:,56) ! wp/dz
1102	hom(:,1,57,sr) = sums(:,57) ! ( w"e + w"p"/rho )/dz
1103	hom(:,1,58,sr) = sums(:,58) ! u"pt"
1104	hom(:,1,59,sr) = sums(:,59) ! upt
1105	hom(:,1,60,sr) = sums(:,58) + sums(:,59) ! upt_t
1106	hom(:,1,61,sr) = sums(:,61) ! v"pt"
1107	hom(:,1,62,sr) = sums(:,62) ! vpt
1108	hom(:,1,63,sr) = sums(:,61) + sums(:,62) ! vpt_t
1109	hom(:,1,64,sr) = sums(:,64) ! rho
1110	hom(:,1,65,sr) = sums(:,65) ! w"sa"
1111	hom(:,1,66,sr) = sums(:,66) ! wsa
1112	hom(:,1,67,sr) = sums(:,65) + sums(:,66) ! wsa
1113	hom(:,1,68,sr) = sums(:,68) ! wp
1114	hom(:,1,69,sr) = sums(:,69) ! w"e + w"p"/rho
1115	hom(:,1,70,sr) = sums(:,70) ! q*2
1116	hom(:,1,71,sr) = sums(:,71) ! prho
1117	hom(:,1,72,sr) = hyp * 1E-4 ! hyp in dbar
1118	hom(:,1,73,sr) = sums(:,73) ! nr
1119	hom(:,1,74,sr) = sums(:,74) ! qr
1120	hom(:,1,75,sr) = sums(:,75) ! qc
1121	hom(:,1,76,sr) = sums(:,76) ! prr (precipitation rate)
1122	! 77 is initial density profile
1123	hom(:,1,78,sr) = ug ! ug
1124	hom(:,1,79,sr) = vg ! vg
1125	hom(:,1,80,sr) = w_subs ! w_subs
1126
1127	hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
1128	! upstream-parts u_x, u_y, u_z, v_x,
1129	! v_y, usw. (in last but one profile)
1130	hom(:,1,pr_palm,sr) = sums(:,pr_palm)
1131	! u, w'u', w'v', t (in last profile)
1132
1133	IF ( max_pr_user > 0 ) THEN ! user-defined profiles
1134	hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
1135	sums(:,pr_palm+1:pr_palm+max_pr_user)
1136	ENDIF
1137
1138	!
1139	!-- Determine the boundary layer height using two different schemes.
1140	!-- First scheme: Starting from the Earth's (Ocean's) surface, look for the
1141	!-- first relative minimum (maximum) of the total heat flux.
1142	!-- The corresponding height is assumed as the boundary layer height, if it
1143	!-- is less than 1.5 times the height where the heat flux becomes negative
1144	!-- (positive) for the first time.
1145	z_i(1) = 0.0
1146	first = .TRUE.
1147
1148	IF ( ocean ) THEN
1149	DO k = nzt, nzb+1, -1
1150	IF ( first .AND. hom(k,1,18,sr) < 0.0 &
1151	.AND. abs(hom(k,1,18,sr)) > 1.0E-8) THEN
1152	first = .FALSE.
1153	height = zw(k)
1154	ENDIF
1155	IF ( hom(k,1,18,sr) < 0.0 .AND. &
1156	abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
1157	hom(k-1,1,18,sr) > hom(k,1,18,sr) ) THEN
1158	IF ( zw(k) < 1.5 * height ) THEN
1159	z_i(1) = zw(k)
1160	ELSE
1161	z_i(1) = height
1162	ENDIF
1163	EXIT
1164	ENDIF
1165	ENDDO
1166	ELSE
1167	DO k = nzb, nzt-1
1168	IF ( first .AND. hom(k,1,18,sr) < 0.0 &
1169	.AND. abs(hom(k,1,18,sr)) > 1.0E-8 ) THEN
1170	first = .FALSE.
1171	height = zw(k)
1172	ENDIF
1173	IF ( hom(k,1,18,sr) < 0.0 .AND. &
1174	abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
1175	hom(k+1,1,18,sr) > hom(k,1,18,sr) ) THEN
1176	IF ( zw(k) < 1.5 * height ) THEN
1177	z_i(1) = zw(k)
1178	ELSE
1179	z_i(1) = height
1180	ENDIF
1181	EXIT
1182	ENDIF
1183	ENDDO
1184	ENDIF
1185
1186	!
1187	!-- Second scheme: Gradient scheme from Sullivan et al. (1998), modified
1188	!-- by Uhlenbrock(2006). The boundary layer height is the height with the
1189	!-- maximal local temperature gradient: starting from the second (the last
1190	!-- but one) vertical gridpoint, the local gradient must be at least
1191	!-- 0.2K/100m and greater than the next four gradients.
1192	!-- WARNING: The threshold value of 0.2K/100m must be adjusted for the
1193	!-- ocean case!
1194	z_i(2) = 0.0
1195	DO k = nzb+1, nzt+1
1196	dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
1197	ENDDO
1198	dptdz_threshold = 0.2_wp / 100.0_wp
1199
1200	IF ( ocean ) THEN
1201	DO k = nzt+1, nzb+5, -1
1202	IF ( dptdz(k) > dptdz_threshold .AND. &
1203	dptdz(k) > dptdz(k-1) .AND. dptdz(k) > dptdz(k-2) .AND. &
1204	dptdz(k) > dptdz(k-3) .AND. dptdz(k) > dptdz(k-4) ) THEN
1205	z_i(2) = zw(k-1)
1206	EXIT
1207	ENDIF
1208	ENDDO
1209	ELSE
1210	DO k = nzb+1, nzt-3
1211	IF ( dptdz(k) > dptdz_threshold .AND. &
1212	dptdz(k) > dptdz(k+1) .AND. dptdz(k) > dptdz(k+2) .AND. &
1213	dptdz(k) > dptdz(k+3) .AND. dptdz(k) > dptdz(k+4) ) THEN
1214	z_i(2) = zw(k-1)
1215	EXIT
1216	ENDIF
1217	ENDDO
1218	ENDIF
1219
1220	hom(nzb+6,1,pr_palm,sr) = z_i(1)
1221	hom(nzb+7,1,pr_palm,sr) = z_i(2)
1222
1223	!
1224	!-- Computation of both the characteristic vertical velocity and
1225	!-- the characteristic convective boundary layer temperature.
1226	!-- The horizontal average at nzb+1 is input for the average temperature.
1227	IF ( hom(nzb,1,18,sr) > 0.0 .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8 &
1228	.AND. z_i(1) /= 0.0 ) THEN
1229	hom(nzb+8,1,pr_palm,sr) = ( g / hom(nzb+1,1,4,sr) * &
1230	hom(nzb,1,18,sr) * &
1231	ABS( z_i(1) ) )**0.333333333
1232	!-- so far this only works if Prandtl layer is used
1233	hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
1234	ELSE
1235	hom(nzb+8,1,pr_palm,sr) = 0.0
1236	hom(nzb+11,1,pr_palm,sr) = 0.0
1237	ENDIF
1238
1239	!
1240	!-- Collect the time series quantities
1241	ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr) ! E
1242	ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr) ! E*
1243	ts_value(3,sr) = dt_3d
1244	ts_value(4,sr) = hom(nzb,1,pr_palm,sr) ! u*
1245	ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr) ! th*
1246	ts_value(6,sr) = u_max
1247	ts_value(7,sr) = v_max
1248	ts_value(8,sr) = w_max
1249	ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr) ! new divergence
1250	ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr) ! old Divergence
1251	ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr) ! z_i(1)
1252	ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr) ! z_i(2)
1253	ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr) ! w*
1254	ts_value(14,sr) = hom(nzb,1,16,sr) ! w'pt' at k=0
1255	ts_value(15,sr) = hom(nzb+1,1,16,sr) ! w'pt' at k=1
1256	ts_value(16,sr) = hom(nzb+1,1,18,sr) ! wpt at k=1
1257	ts_value(17,sr) = hom(nzb,1,4,sr) ! pt(0)
1258	ts_value(18,sr) = hom(nzb+1,1,4,sr) ! pt(zp)
1259	ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr) ! u'w' at k=0
1260	ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr) ! v'w' at k=0
1261	ts_value(21,sr) = hom(nzb,1,48,sr) ! w"q" at k=0
1262
1263	IF ( ts_value(5,sr) /= 0.0 ) THEN
1264	ts_value(22,sr) = ts_value(4,sr)**2 / &
1265	( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
1266	ELSE
1267	ts_value(22,sr) = 10000.0
1268	ENDIF
1269
1270	ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr) ! q*
1271	!
1272	!-- Calculate additional statistics provided by the user interface
1273	CALL user_statistics( 'time_series', sr, 0 )
1274
1275	ENDDO ! loop of the subregions
1276
1277	!
1278	!-- If required, sum up horizontal averages for subsequent time averaging
1279	IF ( do_sum ) THEN
1280	IF ( average_count_pr == 0 ) hom_sum = 0.0
1281	hom_sum = hom_sum + hom(:,1,:,:)
1282	average_count_pr = average_count_pr + 1
1283	do_sum = .FALSE.
1284	ENDIF
1285
1286	!
1287	!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
1288	!-- This flag is reset after each time step in time_integration.
1289	flow_statistics_called = .TRUE.
1290
1291	CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
1292
1293
1294	END SUBROUTINE flow_statistics
1295
1296
1297	#else
1298
1299
1300	!------------------------------------------------------------------------------!
1301	! flow statistics - accelerator version
1302	!------------------------------------------------------------------------------!
1303	SUBROUTINE flow_statistics
1304
1305	USE arrays_3d, &
1306	ONLY: ddzu, ddzw, e, hyp, km, kh,nr, p, prho, pt, q, qc, ql, qr, &
1307	qs, qsws, qswst, rho, sa, saswsb, saswst, shf, ts, tswst, u, &
1308	ug, us, usws, uswst, vsws, v, vg, vpt, vswst, w, w_subs, zw
1309
1310	USE cloud_parameters, &
1311	ONLY: l_d_cp, prr, pt_d_t
1312
1313	USE control_parameters, &
1314	ONLY: average_count_pr, cloud_droplets, cloud_physics, do_sum, &
1315	dt_3d, g, humidity, icloud_scheme, kappa, max_pr_user, &
1316	message_string, ocean, passive_scalar, precipitation, &
1317	use_surface_fluxes, use_top_fluxes, ws_scheme_mom, ws_scheme_sca
1318
1319	USE cpulog, &
1320	ONLY: cpu_log, log_point
1321
1322	USE grid_variables, &
1323	ONLY: ddx, ddy
1324
1325	USE indices, &
1326	ONLY: ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums, nxl, &
1327	nxr, nyn, nys, nzb, nzb_diff_s_inner, nzb_s_inner, nzt, nzt_diff
1328
1329	USE kinds
1330
1331	USE pegrid
1332
1333	USE statistics
1334
1335	IMPLICIT NONE
1336
1337	INTEGER(iwp) :: i !:
1338	INTEGER(iwp) :: j !:
1339	INTEGER(iwp) :: k !:
1340	INTEGER(iwp) :: omp_get_thread_num !:
1341	INTEGER(iwp) :: sr !:
1342	INTEGER(iwp) :: tn !:
1343
1344	LOGICAL :: first !:
1345
1346	REAL(wp) :: dptdz_threshold !:
1347	REAL(wp) :: height !:
1348	REAL(wp) :: pts !:
1349	REAL(wp) :: sums_l_eper !:
1350	REAL(wp) :: sums_l_etot !:
1351	REAL(wp) :: s1 !:
1352	REAL(wp) :: s2 !:
1353	REAL(wp) :: s3 !:
1354	REAL(wp) :: s4 !:
1355	REAL(wp) :: s5 !:
1356	REAL(wp) :: s6 !:
1357	REAL(wp) :: s7 !:
1358	REAL(wp) :: ust !:
1359	REAL(wp) :: ust2 !:
1360	REAL(wp) :: u2, !:
1361	REAL(wp) :: vst !:
1362	REAL(wp) :: vst2 !:
1363	REAL(wp) :: v2 !:
1364	REAL(wp) :: w2 !:
1365	REAL(wp) :: z_i(2) !:
1366
1367	REAL(wp) :: dptdz(nzb+1:nzt+1) !:
1368	REAL(wp) :: sums_ll(nzb:nzt+1,2) !:
1369
1370	CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
1371
1372	!
1373	!-- To be on the safe side, check whether flow_statistics has already been
1374	!-- called once after the current time step
1375	IF ( flow_statistics_called ) THEN
1376
1377	message_string = 'flow_statistics is called two times within one ' // &
1378	'timestep'
1379	CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
1380
1381	ENDIF
1382
1383	!$acc data copyin( hom ) create( sums, sums_l )
1384
1385	!
1386	!-- Compute statistics for each (sub-)region
1387	DO sr = 0, statistic_regions
1388
1389	!
1390	!-- Initialize (local) summation array
1391	sums_l = 0.0
1392
1393	!
1394	!-- Store sums that have been computed in other subroutines in summation
1395	!-- array
1396	sums_l(:,11,:) = sums_l_l(:,sr,:) ! mixing length from diffusivities
1397	!-- WARNING: next line still has to be adjusted for OpenMP
1398	sums_l(:,21,0) = sums_wsts_bc_l(:,sr) ! heat flux from advec_s_bc
1399	sums_l(nzb+9,pr_palm,0) = sums_divold_l(sr) ! old divergence from pres
1400	sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr) ! new divergence from pres
1401
1402	!
1403	!-- Copy the turbulent quantities, evaluated in the advection routines to
1404	!-- the local array sums_l() for further computations
1405	IF ( ws_scheme_mom .AND. sr == 0 ) THEN
1406
1407	!
1408	!-- According to the Neumann bc for the horizontal velocity components,
1409	!-- the corresponding fluxes has to satisfiy the same bc.
1410	IF ( ocean ) THEN
1411	sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
1412	sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)
1413	ENDIF
1414
1415	DO i = 0, threads_per_task-1
1416	!
1417	!-- Swap the turbulent quantities evaluated in advec_ws.
1418	sums_l(:,13,i) = sums_wsus_ws_l(:,i) ! wu
1419	sums_l(:,15,i) = sums_wsvs_ws_l(:,i) ! wv
1420	sums_l(:,30,i) = sums_us2_ws_l(:,i) ! u*2
1421	sums_l(:,31,i) = sums_vs2_ws_l(:,i) ! v*2
1422	sums_l(:,32,i) = sums_ws2_ws_l(:,i) ! w*2
1423	sums_l(:,34,i) = sums_l(:,34,i) + 0.5 * &
1424	( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) + &
1425	sums_ws2_ws_l(:,i) ) ! e*
1426	DO k = nzb, nzt
1427	sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5 * ( &
1428	sums_us2_ws_l(k,i) + &
1429	sums_vs2_ws_l(k,i) + &
1430	sums_ws2_ws_l(k,i) )
1431	ENDDO
1432	ENDDO
1433
1434	ENDIF
1435
1436	IF ( ws_scheme_sca .AND. sr == 0 ) THEN
1437
1438	DO i = 0, threads_per_task-1
1439	sums_l(:,17,i) = sums_wspts_ws_l(:,i) ! wpt from advec_s_ws
1440	IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! wsa
1441	IF ( humidity .OR. passive_scalar ) sums_l(:,49,i) = &
1442	sums_wsqs_ws_l(:,i) !wq
1443	ENDDO
1444
1445	ENDIF
1446	!
1447	!-- Horizontally averaged profiles of horizontal velocities and temperature.
1448	!-- They must have been computed before, because they are already required
1449	!-- for other horizontal averages.
1450	tn = 0
1451
1452	!$OMP PARALLEL PRIVATE( i, j, k, tn )
1453	#if defined( __intel_openmp_bug )
1454	tn = omp_get_thread_num()
1455	#else
1456	!$ tn = omp_get_thread_num()
1457	#endif
1458
1459	!$acc update device( sums_l )
1460
1461	!$OMP DO
1462	!$acc parallel loop gang present( pt, rflags_invers, rmask, sums_l, u, v ) create( s1, s2, s3 )
1463	DO k = nzb, nzt+1
1464	!$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
1465	DO i = nxl, nxr
1466	DO j = nys, nyn
1467	!
1468	!-- k+1 is used in rflags since rflags is set 0 at surface points
1469	s1 = s1 + u(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1470	s2 = s2 + v(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1471	s3 = s3 + pt(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1472	ENDDO
1473	ENDDO
1474	sums_l(k,1,tn) = s1
1475	sums_l(k,2,tn) = s2
1476	sums_l(k,4,tn) = s3
1477	ENDDO
1478	!$acc end parallel loop
1479
1480	!
1481	!-- Horizontally averaged profile of salinity
1482	IF ( ocean ) THEN
1483	!$OMP DO
1484	!$acc parallel loop gang present( rflags_invers, rmask, sums_l, sa ) create( s1 )
1485	DO k = nzb, nzt+1
1486	!$acc loop vector collapse( 2 ) reduction( +: s1 )
1487	DO i = nxl, nxr
1488	DO j = nys, nyn
1489	s1 = s1 + sa(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1490	ENDDO
1491	ENDDO
1492	sums_l(k,23,tn) = s1
1493	ENDDO
1494	!$acc end parallel loop
1495	ENDIF
1496
1497	!
1498	!-- Horizontally averaged profiles of virtual potential temperature,
1499	!-- total water content, specific humidity and liquid water potential
1500	!-- temperature
1501	IF ( humidity ) THEN
1502
1503	!$OMP DO
1504	!$acc parallel loop gang present( q, rflags_invers, rmask, sums_l, vpt ) create( s1, s2 )
1505	DO k = nzb, nzt+1
1506	!$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
1507	DO i = nxl, nxr
1508	DO j = nys, nyn
1509	s1 = s1 + q(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1510	s2 = s2 + vpt(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1511	ENDDO
1512	ENDDO
1513	sums_l(k,41,tn) = s1
1514	sums_l(k,44,tn) = s2
1515	ENDDO
1516	!$acc end parallel loop
1517
1518	IF ( cloud_physics ) THEN
1519	!$OMP DO
1520	!$acc parallel loop gang present( pt, q, ql, rflags_invers, rmask, sums_l ) create( s1, s2 )
1521	DO k = nzb, nzt+1
1522	!$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
1523	DO i = nxl, nxr
1524	DO j = nys, nyn
1525	s1 = s1 + ( q(k,j,i) - ql(k,j,i) ) * &
1526	rmask(j,i,sr) * rflags_invers(j,i,k+1)
1527	s2 = s2 + ( pt(k,j,i) + l_d_cppt_d_t(k) ql(k,j,i) ) * &
1528	rmask(j,i,sr) * rflags_invers(j,i,k+1)
1529	ENDDO
1530	ENDDO
1531	sums_l(k,42,tn) = s1
1532	sums_l(k,43,tn) = s2
1533	ENDDO
1534	!$acc end parallel loop
1535	ENDIF
1536	ENDIF
1537
1538	!
1539	!-- Horizontally averaged profiles of passive scalar
1540	IF ( passive_scalar ) THEN
1541	!$OMP DO
1542	!$acc parallel loop gang present( q, rflags_invers, rmask, sums_l ) create( s1 )
1543	DO k = nzb, nzt+1
1544	!$acc loop vector collapse( 2 ) reduction( +: s1 )
1545	DO i = nxl, nxr
1546	DO j = nys, nyn
1547	s1 = s1 + q(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1548	ENDDO
1549	ENDDO
1550	sums_l(k,41,tn) = s1
1551	ENDDO
1552	!$acc end parallel loop
1553	ENDIF
1554	!$OMP END PARALLEL
1555
1556	!
1557	!-- Summation of thread sums
1558	IF ( threads_per_task > 1 ) THEN
1559	DO i = 1, threads_per_task-1
1560	!$acc parallel present( sums_l )
1561	sums_l(:,1,0) = sums_l(:,1,0) + sums_l(:,1,i)
1562	sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i)
1563	sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i)
1564	!$acc end parallel
1565	IF ( ocean ) THEN
1566	!$acc parallel present( sums_l )
1567	sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)
1568	!$acc end parallel
1569	ENDIF
1570	IF ( humidity ) THEN
1571	!$acc parallel present( sums_l )
1572	sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
1573	sums_l(:,44,0) = sums_l(:,44,0) + sums_l(:,44,i)
1574	!$acc end parallel
1575	IF ( cloud_physics ) THEN
1576	!$acc parallel present( sums_l )
1577	sums_l(:,42,0) = sums_l(:,42,0) + sums_l(:,42,i)
1578	sums_l(:,43,0) = sums_l(:,43,0) + sums_l(:,43,i)
1579	!$acc end parallel
1580	ENDIF
1581	ENDIF
1582	IF ( passive_scalar ) THEN
1583	!$acc parallel present( sums_l )
1584	sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
1585	!$acc end parallel
1586	ENDIF
1587	ENDDO
1588	ENDIF
1589
1590	#if defined( __parallel )
1591	!
1592	!-- Compute total sum from local sums
1593	!$acc update host( sums_l )
1594	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1595	CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL, &
1596	MPI_SUM, comm2d, ierr )
1597	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1598	CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL, &
1599	MPI_SUM, comm2d, ierr )
1600	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1601	CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL, &
1602	MPI_SUM, comm2d, ierr )
1603	IF ( ocean ) THEN
1604	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1605	CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb, &
1606	MPI_REAL, MPI_SUM, comm2d, ierr )
1607	ENDIF
1608	IF ( humidity ) THEN
1609	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1610	CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb, &
1611	MPI_REAL, MPI_SUM, comm2d, ierr )
1612	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1613	CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
1614	MPI_REAL, MPI_SUM, comm2d, ierr )
1615	IF ( cloud_physics ) THEN
1616	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1617	CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb, &
1618	MPI_REAL, MPI_SUM, comm2d, ierr )
1619	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1620	CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb, &
1621	MPI_REAL, MPI_SUM, comm2d, ierr )
1622	ENDIF
1623	ENDIF
1624
1625	IF ( passive_scalar ) THEN
1626	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1627	CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
1628	MPI_REAL, MPI_SUM, comm2d, ierr )
1629	ENDIF
1630	!$acc update device( sums )
1631	#else
1632	!$acc parallel present( sums, sums_l )
1633	sums(:,1) = sums_l(:,1,0)
1634	sums(:,2) = sums_l(:,2,0)
1635	sums(:,4) = sums_l(:,4,0)
1636	!$acc end parallel
1637	IF ( ocean ) THEN
1638	!$acc parallel present( sums, sums_l )
1639	sums(:,23) = sums_l(:,23,0)
1640	!$acc end parallel
1641	ENDIF
1642	IF ( humidity ) THEN
1643	!$acc parallel present( sums, sums_l )
1644	sums(:,44) = sums_l(:,44,0)
1645	sums(:,41) = sums_l(:,41,0)
1646	!$acc end parallel
1647	IF ( cloud_physics ) THEN
1648	!$acc parallel present( sums, sums_l )
1649	sums(:,42) = sums_l(:,42,0)
1650	sums(:,43) = sums_l(:,43,0)
1651	!$acc end parallel
1652	ENDIF
1653	ENDIF
1654	IF ( passive_scalar ) THEN
1655	!$acc parallel present( sums, sums_l )
1656	sums(:,41) = sums_l(:,41,0)
1657	!$acc end parallel
1658	ENDIF
1659	#endif
1660
1661	!
1662	!-- Final values are obtained by division by the total number of grid points
1663	!-- used for summation. After that store profiles.
1664	!$acc parallel present( hom, ngp_2dh, ngp_2dh_s_inner, sums )
1665	sums(:,1) = sums(:,1) / ngp_2dh(sr)
1666	sums(:,2) = sums(:,2) / ngp_2dh(sr)
1667	sums(:,4) = sums(:,4) / ngp_2dh_s_inner(:,sr)
1668	hom(:,1,1,sr) = sums(:,1) ! u
1669	hom(:,1,2,sr) = sums(:,2) ! v
1670	hom(:,1,4,sr) = sums(:,4) ! pt
1671	!$acc end parallel
1672
1673	!
1674	!-- Salinity
1675	IF ( ocean ) THEN
1676	!$acc parallel present( hom, ngp_2dh_s_inner, sums )
1677	sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr)
1678	hom(:,1,23,sr) = sums(:,23) ! sa
1679	!$acc end parallel
1680	ENDIF
1681
1682	!
1683	!-- Humidity and cloud parameters
1684	IF ( humidity ) THEN
1685	!$acc parallel present( hom, ngp_2dh_s_inner, sums )
1686	sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr)
1687	sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
1688	hom(:,1,44,sr) = sums(:,44) ! vpt
1689	hom(:,1,41,sr) = sums(:,41) ! qv (q)
1690	!$acc end parallel
1691	IF ( cloud_physics ) THEN
1692	!$acc parallel present( hom, ngp_2dh_s_inner, sums )
1693	sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr)
1694	sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr)
1695	hom(:,1,42,sr) = sums(:,42) ! qv
1696	hom(:,1,43,sr) = sums(:,43) ! pt
1697	!$acc end parallel
1698	ENDIF
1699	ENDIF
1700
1701	!
1702	!-- Passive scalar
1703	IF ( passive_scalar ) THEN
1704	!$acc parallel present( hom, ngp_2dh_s_inner, sums )
1705	sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
1706	hom(:,1,41,sr) = sums(:,41) ! s (q)
1707	!$acc end parallel
1708	ENDIF
1709
1710	!
1711	!-- Horizontally averaged profiles of the remaining prognostic variables,
1712	!-- variances, the total and the perturbation energy (single values in last
1713	!-- column of sums_l) and some diagnostic quantities.
1714	!-- NOTE: for simplicity, nzb_s_inner is used below, although strictly
1715	!-- ---- speaking the following k-loop would have to be split up and
1716	!-- rearranged according to the staggered grid.
1717	!-- However, this implies no error since staggered velocity components
1718	!-- are zero at the walls and inside buildings.
1719	tn = 0
1720	#if defined( __intel_openmp_bug )
1721	!$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, &
1722	!$OMP tn, ust, ust2, u2, vst, vst2, v2, w2 )
1723	tn = omp_get_thread_num()
1724	#else
1725	!$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2, w2 )
1726	!$ tn = omp_get_thread_num()
1727	#endif
1728	!$OMP DO
1729	!$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 )
1730	DO k = nzb, nzt+1
1731	!$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5, s6, s7 )
1732	DO i = nxl, nxr
1733	DO j = nys, nyn
1734	!
1735	!-- Prognostic and diagnostic variables
1736	s1 = s1 + w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1737	s2 = s2 + e(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1738	s3 = s3 + km(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1739	s4 = s4 + kh(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1740	s5 = s5 + p(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1741	s6 = s6 + ( pt(k,j,i)-hom(k,1,4,sr) )*2 rmask(j,i,sr) * &
1742	rflags_invers(j,i,k+1)
1743	!
1744	!-- Higher moments
1745	!-- (Computation of the skewness of w further below)
1746	s7 = s7 + w(k,j,i)*3 rmask(j,i,sr) * rflags_invers(j,i,k+1)
1747	ENDDO
1748	ENDDO
1749	sums_l(k,3,tn) = s1
1750	sums_l(k,8,tn) = s2
1751	sums_l(k,9,tn) = s3
1752	sums_l(k,10,tn) = s4
1753	sums_l(k,40,tn) = s5
1754	sums_l(k,33,tn) = s6
1755	sums_l(k,38,tn) = s7
1756	ENDDO
1757	!$acc end parallel loop
1758
1759	IF ( humidity ) THEN
1760	!$OMP DO
1761	!$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l ) create( s1 )
1762	DO k = nzb, nzt+1
1763	!$acc loop vector collapse( 2 ) reduction( +: s1 )
1764	DO i = nxl, nxr
1765	DO j = nys, nyn
1766	s1 = s1 + ( q(k,j,i)-hom(k,1,41,sr) )*2 rmask(j,i,sr) * &
1767	rflags_invers(j,i,k+1)
1768	ENDDO
1769	ENDDO
1770	sums_l(k,70,tn) = s1
1771	ENDDO
1772	!$acc end parallel loop
1773	ENDIF
1774
1775	!
1776	!-- Total and perturbation energy for the total domain (being
1777	!-- collected in the last column of sums_l).
1778	!$OMP DO
1779	!$acc parallel loop collapse(3) present( rflags_invers, rmask, u, v, w ) reduction(+:s1)
1780	DO i = nxl, nxr
1781	DO j = nys, nyn
1782	DO k = nzb, nzt+1
1783	s1 = s1 + 0.5 * ( u(k,j,i)2 + v(k,j,i)2 + w(k,j,i)*2 ) &
1784	rmask(j,i,sr) * rflags_invers(j,i,k+1)
1785	ENDDO
1786	ENDDO
1787	ENDDO
1788	!$acc end parallel loop
1789	!$acc parallel present( sums_l )
1790	sums_l(nzb+4,pr_palm,tn) = s1
1791	!$acc end parallel
1792
1793	!$OMP DO
1794	!$acc parallel present( rmask, sums_l, us, usws, vsws, ts ) create( s1, s2, s3, s4 )
1795	!$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4 )
1796	DO i = nxl, nxr
1797	DO j = nys, nyn
1798	!
1799	!-- 2D-arrays (being collected in the last column of sums_l)
1800	s1 = s1 + us(j,i) * rmask(j,i,sr)
1801	s2 = s2 + usws(j,i) * rmask(j,i,sr)
1802	s3 = s3 + vsws(j,i) * rmask(j,i,sr)
1803	s4 = s4 + ts(j,i) * rmask(j,i,sr)
1804	ENDDO
1805	ENDDO
1806	sums_l(nzb,pr_palm,tn) = s1
1807	sums_l(nzb+1,pr_palm,tn) = s2
1808	sums_l(nzb+2,pr_palm,tn) = s3
1809	sums_l(nzb+3,pr_palm,tn) = s4
1810	!$acc end parallel
1811
1812	IF ( humidity ) THEN
1813	!$acc parallel present( qs, rmask, sums_l ) create( s1 )
1814	!$acc loop vector collapse( 2 ) reduction( +: s1 )
1815	DO i = nxl, nxr
1816	DO j = nys, nyn
1817	s1 = s1 + qs(j,i) * rmask(j,i,sr)
1818	ENDDO
1819	ENDDO
1820	sums_l(nzb+12,pr_palm,tn) = s1
1821	!$acc end parallel
1822	ENDIF
1823
1824	!
1825	!-- Computation of statistics when ws-scheme is not used. Else these
1826	!-- quantities are evaluated in the advection routines.
1827	IF ( .NOT. ws_scheme_mom .OR. sr /= 0 ) THEN
1828
1829	!$OMP DO
1830	!$acc parallel loop gang present( u, v, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3, s4, ust2, vst2, w2 )
1831	DO k = nzb, nzt+1
1832	!$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4 )
1833	DO i = nxl, nxr
1834	DO j = nys, nyn
1835	ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
1836	vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
1837	w2 = w(k,j,i)**2
1838
1839	s1 = s1 + ust2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1840	s2 = s2 + vst2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1841	s3 = s3 + w2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1842	!
1843	!-- Perturbation energy
1844	s4 = s4 + 0.5 * ( ust2 + vst2 + w2 ) * rmask(j,i,sr) * &
1845	rflags_invers(j,i,k+1)
1846	ENDDO
1847	ENDDO
1848	sums_l(k,30,tn) = s1
1849	sums_l(k,31,tn) = s2
1850	sums_l(k,32,tn) = s3
1851	sums_l(k,34,tn) = s4
1852	ENDDO
1853	!$acc end parallel loop
1854	!
1855	!-- Total perturbation TKE
1856	!$OMP DO
1857	!$acc parallel present( sums_l ) create( s1 )
1858	!$acc loop reduction( +: s1 )
1859	DO k = nzb, nzt+1
1860	s1 = s1 + sums_l(k,34,tn)
1861	ENDDO
1862	sums_l(nzb+5,pr_palm,tn) = s1
1863	!$acc end parallel
1864
1865	ENDIF
1866
1867	!
1868	!-- Horizontally averaged profiles of the vertical fluxes
1869
1870	!
1871	!-- Subgridscale fluxes.
1872	!-- WARNING: If a Prandtl-layer is used (k=nzb for flat terrain), the fluxes
1873	!-- ------- should be calculated there in a different way. This is done
1874	!-- in the next loop further below, where results from this loop are
1875	!-- overwritten. However, THIS WORKS IN CASE OF FLAT TERRAIN ONLY!
1876	!-- The non-flat case still has to be handled.
1877	!-- NOTE: for simplicity, nzb_s_inner is used below, although
1878	!-- ---- strictly speaking the following k-loop would have to be
1879	!-- split up according to the staggered grid.
1880	!-- However, this implies no error since staggered velocity
1881	!-- components are zero at the walls and inside buildings.
1882	!$OMP DO
1883	!$acc parallel loop gang present( ddzu, kh, km, pt, u, v, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3 )
1884	DO k = nzb, nzt_diff
1885	!$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
1886	DO i = nxl, nxr
1887	DO j = nys, nyn
1888
1889	!
1890	!-- Momentum flux w"u"
1891	s1 = s1 - 0.25 * ( &
1892	km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
1893	) * ( &
1894	( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
1895	+ ( w(k,j,i) - w(k,j,i-1) ) * ddx &
1896	) &
1897	* rmask(j,i,sr) * rflags_invers(j,i,k+1)
1898	!
1899	!-- Momentum flux w"v"
1900	s2 = s2 - 0.25 * ( &
1901	km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
1902	) * ( &
1903	( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
1904	+ ( w(k,j,i) - w(k,j-1,i) ) * ddy &
1905	) &
1906	* rmask(j,i,sr) * rflags_invers(j,i,k+1)
1907	!
1908	!-- Heat flux w"pt"
1909	s3 = s3 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) ) &
1910	* ( pt(k+1,j,i) - pt(k,j,i) ) &
1911	* ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1912	ENDDO
1913	ENDDO
1914	sums_l(k,12,tn) = s1
1915	sums_l(k,14,tn) = s2
1916	sums_l(k,16,tn) = s3
1917	ENDDO
1918	!$acc end parallel loop
1919
1920	!
1921	!-- Salinity flux w"sa"
1922	IF ( ocean ) THEN
1923	!$acc parallel loop gang present( ddzu, kh, sa, rflags_invers, rmask, sums_l ) create( s1 )
1924	DO k = nzb, nzt_diff
1925	!$acc loop vector collapse( 2 ) reduction( +: s1 )
1926	DO i = nxl, nxr
1927	DO j = nys, nyn
1928	s1 = s1 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) ) &
1929	* ( sa(k+1,j,i) - sa(k,j,i) ) &
1930	* ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1931	ENDDO
1932	ENDDO
1933	sums_l(k,65,tn) = s1
1934	ENDDO
1935	!$acc end parallel loop
1936	ENDIF
1937
1938	!
1939	!-- Buoyancy flux, water flux (humidity flux) w"q"
1940	IF ( humidity ) THEN
1941
1942	!$acc parallel loop gang present( ddzu, kh, q, vpt, rflags_invers, rmask, sums_l ) create( s1, s2 )
1943	DO k = nzb, nzt_diff
1944	!$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
1945	DO i = nxl, nxr
1946	DO j = nys, nyn
1947	s1 = s1 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) ) &
1948	* ( vpt(k+1,j,i) - vpt(k,j,i) ) &
1949	* ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1950	s2 = s2 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) ) &
1951	* ( q(k+1,j,i) - q(k,j,i) ) &
1952	* ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1953	ENDDO
1954	ENDDO
1955	sums_l(k,45,tn) = s1
1956	sums_l(k,48,tn) = s2
1957	ENDDO
1958	!$acc end parallel loop
1959
1960	IF ( cloud_physics ) THEN
1961
1962	!$acc parallel loop gang present( ddzu, kh, q, ql, rflags_invers, rmask, sums_l ) create( s1 )
1963	DO k = nzb, nzt_diff
1964	!$acc loop vector collapse( 2 ) reduction( +: s1 )
1965	DO i = nxl, nxr
1966	DO j = nys, nyn
1967	s1 = s1 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) ) &
1968	* ( ( q(k+1,j,i) - ql(k+1,j,i) ) &
1969	- ( q(k,j,i) - ql(k,j,i) ) ) &
1970	* ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1971	ENDDO
1972	ENDDO
1973	sums_l(k,51,tn) = s1
1974	ENDDO
1975	!$acc end parallel loop
1976
1977	ENDIF
1978
1979	ENDIF
1980	!
1981	!-- Passive scalar flux
1982	IF ( passive_scalar ) THEN
1983
1984	!$acc parallel loop gang present( ddzu, kh, q, rflags_invers, rmask, sums_l ) create( s1 )
1985	DO k = nzb, nzt_diff
1986	!$acc loop vector collapse( 2 ) reduction( +: s1 )
1987	DO i = nxl, nxr
1988	DO j = nys, nyn
1989	s1 = s1 - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) ) &
1990	* ( q(k+1,j,i) - q(k,j,i) ) &
1991	* ddzu(k+1) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1992	ENDDO
1993	ENDDO
1994	sums_l(k,48,tn) = s1
1995	ENDDO
1996	!$acc end parallel loop
1997
1998	ENDIF
1999
2000	IF ( use_surface_fluxes ) THEN
2001
2002	!$OMP DO
2003	!$acc parallel present( rmask, shf, sums_l, usws, vsws ) create( s1, s2, s3, s4, s5 )
2004	!$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5 )
2005	DO i = nxl, nxr
2006	DO j = nys, nyn
2007	!
2008	!-- Subgridscale fluxes in the Prandtl layer
2009	s1 = s1 + usws(j,i) * rmask(j,i,sr) ! w"u"
2010	s2 = s2 + vsws(j,i) * rmask(j,i,sr) ! w"v"
2011	s3 = s3 + shf(j,i) * rmask(j,i,sr) ! w"pt"
2012	s4 = s4 + 0.0 * rmask(j,i,sr) ! u"pt"
2013	s5 = s5 + 0.0 * rmask(j,i,sr) ! v"pt"
2014	ENDDO
2015	ENDDO
2016	sums_l(nzb,12,tn) = s1
2017	sums_l(nzb,14,tn) = s2
2018	sums_l(nzb,16,tn) = s3
2019	sums_l(nzb,58,tn) = s4
2020	sums_l(nzb,61,tn) = s5
2021	!$acc end parallel
2022
2023	IF ( ocean ) THEN
2024
2025	!$OMP DO
2026	!$acc parallel present( rmask, saswsb, sums_l ) create( s1 )
2027	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2028	DO i = nxl, nxr
2029	DO j = nys, nyn
2030	s1 = s1 + saswsb(j,i) * rmask(j,i,sr) ! w"sa"
2031	ENDDO
2032	ENDDO
2033	sums_l(nzb,65,tn) = s1
2034	!$acc end parallel
2035
2036	ENDIF
2037
2038	IF ( humidity ) THEN
2039
2040	!$OMP DO
2041	!$acc parallel present( pt, q, qsws, rmask, shf, sums_l ) create( s1, s2 )
2042	!$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2043	DO i = nxl, nxr
2044	DO j = nys, nyn
2045	s1 = s1 + qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
2046	s2 = s2 + ( ( 1.0 + 0.61 * q(nzb,j,i) ) * shf(j,i) &
2047	+ 0.61 * pt(nzb,j,i) * qsws(j,i) )
2048	ENDDO
2049	ENDDO
2050	sums_l(nzb,48,tn) = s1
2051	sums_l(nzb,45,tn) = s2
2052	!$acc end parallel
2053
2054	IF ( cloud_droplets ) THEN
2055
2056	!$OMP DO
2057	!$acc parallel present( pt, q, ql, qsws, rmask, shf, sums_l ) create( s1 )
2058	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2059	DO i = nxl, nxr
2060	DO j = nys, nyn
2061	s1 = s1 + ( ( 1.0 + 0.61 * q(nzb,j,i) - ql(nzb,j,i) ) * &
2062	shf(j,i) + 0.61 * pt(nzb,j,i) * qsws(j,i) )
2063	ENDDO
2064	ENDDO
2065	sums_l(nzb,45,tn) = s1
2066	!$acc end parallel
2067
2068	ENDIF
2069
2070	IF ( cloud_physics ) THEN
2071
2072	!$OMP DO
2073	!$acc parallel present( qsws, rmask, sums_l ) create( s1 )
2074	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2075	DO i = nxl, nxr
2076	DO j = nys, nyn
2077	!
2078	!-- Formula does not work if ql(nzb) /= 0.0
2079	s1 = s1 + qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
2080	ENDDO
2081	ENDDO
2082	sums_l(nzb,51,tn) = s1
2083	!$acc end parallel
2084
2085	ENDIF
2086
2087	ENDIF
2088
2089	IF ( passive_scalar ) THEN
2090
2091	!$OMP DO
2092	!$acc parallel present( qsws, rmask, sums_l ) create( s1 )
2093	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2094	DO i = nxl, nxr
2095	DO j = nys, nyn
2096	s1 = s1 + qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
2097	ENDDO
2098	ENDDO
2099	sums_l(nzb,48,tn) = s1
2100	!$acc end parallel
2101
2102	ENDIF
2103
2104	ENDIF
2105
2106	!
2107	!-- Subgridscale fluxes at the top surface
2108	IF ( use_top_fluxes ) THEN
2109
2110	!$OMP DO
2111	!$acc parallel present( rmask, sums_l, tswst, uswst, vswst ) create( s1, s2, s3, s4, s5 )
2112	!$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5 )
2113	DO i = nxl, nxr
2114	DO j = nys, nyn
2115	s1 = s1 + uswst(j,i) * rmask(j,i,sr) ! w"u"
2116	s2 = s2 + vswst(j,i) * rmask(j,i,sr) ! w"v"
2117	s3 = s3 + tswst(j,i) * rmask(j,i,sr) ! w"pt"
2118	s4 = s4 + 0.0 * rmask(j,i,sr) ! u"pt"
2119	s5 = s5 + 0.0 * rmask(j,i,sr) ! v"pt"
2120	ENDDO
2121	ENDDO
2122	sums_l(nzt:nzt+1,12,tn) = s1
2123	sums_l(nzt:nzt+1,14,tn) = s2
2124	sums_l(nzt:nzt+1,16,tn) = s3
2125	sums_l(nzt:nzt+1,58,tn) = s4
2126	sums_l(nzt:nzt+1,61,tn) = s5
2127	!$acc end parallel
2128
2129	IF ( ocean ) THEN
2130
2131	!$OMP DO
2132	!$acc parallel present( rmask, saswst, sums_l ) create( s1 )
2133	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2134	DO i = nxl, nxr
2135	DO j = nys, nyn
2136	s1 = s1 + saswst(j,i) * rmask(j,i,sr) ! w"sa"
2137	ENDDO
2138	ENDDO
2139	sums_l(nzt,65,tn) = s1
2140	!$acc end parallel
2141
2142	ENDIF
2143
2144	IF ( humidity ) THEN
2145
2146	!$OMP DO
2147	!$acc parallel present( pt, q, qswst, rmask, tswst, sums_l ) create( s1, s2 )
2148	!$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2149	DO i = nxl, nxr
2150	DO j = nys, nyn
2151	s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
2152	s2 = s2 + ( ( 1.0 + 0.61 * q(nzt,j,i) ) * tswst(j,i) + &
2153	0.61 * pt(nzt,j,i) * qswst(j,i) )
2154	ENDDO
2155	ENDDO
2156	sums_l(nzt,48,tn) = s1
2157	sums_l(nzt,45,tn) = s2
2158	!$acc end parallel
2159
2160	IF ( cloud_droplets ) THEN
2161
2162	!$OMP DO
2163	!$acc parallel present( pt, q, ql, qswst, rmask, tswst, sums_l ) create( s1 )
2164	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2165	DO i = nxl, nxr
2166	DO j = nys, nyn
2167	s1 = s1 + ( ( 1.0 + 0.61 * q(nzt,j,i) - ql(nzt,j,i) ) * &
2168	tswst(j,i) + 0.61 * pt(nzt,j,i) * qswst(j,i) )
2169	ENDDO
2170	ENDDO
2171	sums_l(nzt,45,tn) = s1
2172	!$acc end parallel
2173
2174	ENDIF
2175
2176	IF ( cloud_physics ) THEN
2177
2178	!$OMP DO
2179	!$acc parallel present( qswst, rmask, sums_l ) create( s1 )
2180	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2181	DO i = nxl, nxr
2182	DO j = nys, nyn
2183	!
2184	!-- Formula does not work if ql(nzb) /= 0.0
2185	s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
2186	ENDDO
2187	ENDDO
2188	sums_l(nzt,51,tn) = s1
2189	!$acc end parallel
2190
2191	ENDIF
2192
2193	ENDIF
2194
2195	IF ( passive_scalar ) THEN
2196
2197	!$OMP DO
2198	!$acc parallel present( qswst, rmask, sums_l ) create( s1 )
2199	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2200	DO i = nxl, nxr
2201	DO j = nys, nyn
2202	s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
2203	ENDDO
2204	ENDDO
2205	sums_l(nzt,48,tn) = s1
2206	!$acc end parallel
2207
2208	ENDIF
2209
2210	ENDIF
2211
2212	!
2213	!-- Resolved fluxes (can be computed for all horizontal points)
2214	!-- NOTE: for simplicity, nzb_s_inner is used below, although strictly
2215	!-- ---- speaking the following k-loop would have to be split up and
2216	!-- rearranged according to the staggered grid.
2217	!$acc parallel loop gang present( hom, pt, rflags_invers, rmask, sums_l, u, v, w ) create( s1, s2, s3 )
2218	DO k = nzb, nzt_diff
2219	!$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
2220	DO i = nxl, nxr
2221	DO j = nys, nyn
2222	ust = 0.5 * ( u(k,j,i) - hom(k,1,1,sr) + &
2223	u(k+1,j,i) - hom(k+1,1,1,sr) )
2224	vst = 0.5 * ( v(k,j,i) - hom(k,1,2,sr) + &
2225	v(k+1,j,i) - hom(k+1,1,2,sr) )
2226	pts = 0.5 * ( pt(k,j,i) - hom(k,1,4,sr) + &
2227	pt(k+1,j,i) - hom(k+1,1,4,sr) )
2228	!
2229	!-- Higher moments
2230	s1 = s1 + pts * w(k,j,i)*2 rmask(j,i,sr) * rflags_invers(j,i,k+1)
2231	s2 = s2 + pts*2 w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2232	!
2233	!-- Energy flux we (has to be adjusted?)
2234	s3 = s3 + w(k,j,i) * 0.5 * ( ust2 + vst2 + w(k,j,i)**2 ) &
2235	* rmask(j,i,sr) * rflags_invers(j,i,k+1)
2236	ENDDO
2237	ENDDO
2238	sums_l(k,35,tn) = s1
2239	sums_l(k,36,tn) = s2
2240	sums_l(k,37,tn) = s3
2241	ENDDO
2242	!$acc end parallel loop
2243
2244	!
2245	!-- Salinity flux and density (density does not belong to here,
2246	!-- but so far there is no other suitable place to calculate)
2247	IF ( ocean ) THEN
2248
2249	IF( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN
2250
2251	!$acc parallel loop gang present( hom, rflags_invers, rmask, sa, sums_l, w ) create( s1 )
2252	DO k = nzb, nzt_diff
2253	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2254	DO i = nxl, nxr
2255	DO j = nys, nyn
2256	s1 = s1 + 0.5 * ( sa(k,j,i) - hom(k,1,23,sr) + &
2257	sa(k+1,j,i) - hom(k+1,1,23,sr) ) &
2258	* w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2259	ENDDO
2260	ENDDO
2261	sums_l(k,66,tn) = s1
2262	ENDDO
2263	!$acc end parallel loop
2264
2265	ENDIF
2266
2267	!$acc parallel loop gang present( rflags_invers, rho, prho, rmask, sums_l ) create( s1, s2 )
2268	DO k = nzb, nzt_diff
2269	!$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2270	DO i = nxl, nxr
2271	DO j = nys, nyn
2272	s1 = s1 + rho(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2273	s2 = s2 + prho(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2274	ENDDO
2275	ENDDO
2276	sums_l(k,64,tn) = s1
2277	sums_l(k,71,tn) = s2
2278	ENDDO
2279	!$acc end parallel loop
2280
2281	ENDIF
2282
2283	!
2284	!-- Buoyancy flux, water flux, humidity flux, liquid water
2285	!-- content, rain drop concentration and rain water content
2286	IF ( humidity ) THEN
2287
2288	IF ( cloud_physics .OR. cloud_droplets ) THEN
2289
2290	!$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, vpt, w ) create( s1 )
2291	DO k = nzb, nzt_diff
2292	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2293	DO i = nxl, nxr
2294	DO j = nys, nyn
2295	s1 = s1 + 0.5 * ( vpt(k,j,i) - hom(k,1,44,sr) + &
2296	vpt(k+1,j,i) - hom(k+1,1,44,sr) ) * &
2297	w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2298	ENDDO
2299	ENDDO
2300	sums_l(k,46,tn) = s1
2301	ENDDO
2302	!$acc end parallel loop
2303
2304	IF ( .NOT. cloud_droplets ) THEN
2305
2306	!$acc parallel loop gang present( hom, q, ql, rflags_invers, rmask, sums_l, w ) create( s1 )
2307	DO k = nzb, nzt_diff
2308	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2309	DO i = nxl, nxr
2310	DO j = nys, nyn
2311	s1 = s1 + 0.5 * ( ( q(k,j,i) - ql(k,j,i) ) - hom(k,1,42,sr) + &
2312	( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) ) &
2313	* w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2314	ENDDO
2315	ENDDO
2316	sums_l(k,52,tn) = s1
2317	ENDDO
2318	!$acc end parallel loop
2319
2320	IF ( icloud_scheme == 0 ) THEN
2321
2322	!$acc parallel loop gang present( qc, ql, rflags_invers, rmask, sums_l ) create( s1, s2 )
2323	DO k = nzb, nzt_diff
2324	!$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2325	DO i = nxl, nxr
2326	DO j = nys, nyn
2327	s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2328	s2 = s2 + qc(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2329	ENDDO
2330	ENDDO
2331	sums_l(k,54,tn) = s1
2332	sums_l(k,75,tn) = s2
2333	ENDDO
2334	!$acc end parallel loop
2335
2336	IF ( precipitation ) THEN
2337
2338	!$acc parallel loop gang present( nr, qr, prr, rflags_invers, rmask, sums_l ) create( s1, s2, s3 )
2339	DO k = nzb, nzt_diff
2340	!$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
2341	DO i = nxl, nxr
2342	DO j = nys, nyn
2343	s1 = s1 + nr(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2344	s2 = s2 + qr(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2345	s3 = s3 + prr(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2346	ENDDO
2347	ENDDO
2348	sums_l(k,73,tn) = s1
2349	sums_l(k,74,tn) = s2
2350	sums_l(k,76,tn) = s3
2351	ENDDO
2352	!$acc end parallel loop
2353
2354	ENDIF
2355
2356	ELSE
2357
2358	!$acc parallel loop gang present( ql, rflags_invers, rmask, sums_l ) create( s1 )
2359	DO k = nzb, nzt_diff
2360	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2361	DO i = nxl, nxr
2362	DO j = nys, nyn
2363	s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2364	ENDDO
2365	ENDDO
2366	sums_l(k,54,tn) = s1
2367	ENDDO
2368	!$acc end parallel loop
2369
2370	ENDIF
2371
2372	ELSE
2373
2374	!$acc parallel loop gang present( ql, rflags_invers, rmask, sums_l ) create( s1 )
2375	DO k = nzb, nzt_diff
2376	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2377	DO i = nxl, nxr
2378	DO j = nys, nyn
2379	s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2380	ENDDO
2381	ENDDO
2382	sums_l(k,54,tn) = s1
2383	ENDDO
2384	!$acc end parallel loop
2385
2386	ENDIF
2387
2388	ELSE
2389
2390	IF( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN
2391
2392	!$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, vpt, w ) create( s1 )
2393	DO k = nzb, nzt_diff
2394	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2395	DO i = nxl, nxr
2396	DO j = nys, nyn
2397	s1 = s1 + 0.5 * ( vpt(k,j,i) - hom(k,1,44,sr) + &
2398	vpt(k+1,j,i) - hom(k+1,1,44,sr) ) &
2399	* w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2400	ENDDO
2401	ENDDO
2402	sums_l(k,46,tn) = s1
2403	ENDDO
2404	!$acc end parallel loop
2405
2406	ELSEIF ( ws_scheme_sca .AND. sr == 0 ) THEN
2407
2408	!$acc parallel loop present( hom, sums_l )
2409	DO k = nzb, nzt_diff
2410	sums_l(k,46,tn) = ( 1.0 + 0.61 * hom(k,1,41,sr) ) * sums_l(k,17,tn) + &
2411	0.61 * hom(k,1,4,sr) * sums_l(k,49,tn)
2412	ENDDO
2413	!$acc end parallel loop
2414
2415	ENDIF
2416
2417	ENDIF
2418
2419	ENDIF
2420	!
2421	!-- Passive scalar flux
2422	IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca .OR. sr /= 0 ) ) THEN
2423
2424	!$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l, w ) create( s1 )
2425	DO k = nzb, nzt_diff
2426	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2427	DO i = nxl, nxr
2428	DO j = nys, nyn
2429	s1 = s1 + 0.5 * ( q(k,j,i) - hom(k,1,41,sr) + &
2430	q(k+1,j,i) - hom(k+1,1,41,sr) ) &
2431	* w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2432	ENDDO
2433	ENDDO
2434	sums_l(k,49,tn) = s1
2435	ENDDO
2436	!$acc end parallel loop
2437
2438	ENDIF
2439
2440	!
2441	!-- For speed optimization fluxes which have been computed in part directly
2442	!-- inside the WS advection routines are treated seperatly
2443	!-- Momentum fluxes first:
2444	IF ( .NOT. ws_scheme_mom .OR. sr /= 0 ) THEN
2445
2446	!$OMP DO
2447	!$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, u, v, w ) create( s1, s2 )
2448	DO k = nzb, nzt_diff
2449	!$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2450	DO i = nxl, nxr
2451	DO j = nys, nyn
2452	ust = 0.5 * ( u(k,j,i) - hom(k,1,1,sr) + &
2453	u(k+1,j,i) - hom(k+1,1,1,sr) )
2454	vst = 0.5 * ( v(k,j,i) - hom(k,1,2,sr) + &
2455	v(k+1,j,i) - hom(k+1,1,2,sr) )
2456	!
2457	!-- Momentum flux wu
2458	s1 = s1 + 0.5 * ( w(k,j,i-1) + w(k,j,i) ) &
2459	* ust * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2460	!
2461	!-- Momentum flux wv
2462	s2 = s2 + 0.5 * ( w(k,j-1,i) + w(k,j,i) ) &
2463	* vst * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2464	ENDDO
2465	ENDDO
2466	sums_l(k,13,tn) = s1
2467	sums_l(k,15,tn) = s1
2468	ENDDO
2469	!$acc end parallel loop
2470
2471	ENDIF
2472
2473	IF ( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN
2474
2475	!$OMP DO
2476	!$acc parallel loop gang present( hom, pt, rflags_invers, rmask, sums_l, w ) create( s1 )
2477	DO k = nzb, nzt_diff
2478	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2479	DO i = nxl, nxr
2480	DO j = nys, nyn
2481	!
2482	!-- Vertical heat flux
2483	s1 = s1 + 0.5 * ( pt(k,j,i) - hom(k,1,4,sr) + &
2484	pt(k+1,j,i) - hom(k+1,1,4,sr) ) &
2485	* w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2486	ENDDO
2487	ENDDO
2488	sums_l(k,17,tn) = s1
2489	ENDDO
2490	!$acc end parallel loop
2491
2492	IF ( humidity ) THEN
2493
2494	!$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l, w ) create( s1 )
2495	DO k = nzb, nzt_diff
2496	!$acc loop vector collapse( 2 ) reduction( +: s1 )
2497	DO i = nxl, nxr
2498	DO j = nys, nyn
2499	s1 = s1 + 0.5 * ( q(k,j,i) - hom(k,1,41,sr) + &
2500	q(k+1,j,i) - hom(k+1,1,41,sr) ) &
2501	* w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2502	ENDDO
2503	ENDDO
2504	sums_l(k,49,tn) = s1
2505	ENDDO
2506	!$acc end parallel loop
2507
2508	ENDIF
2509
2510	ENDIF
2511
2512
2513	!
2514	!-- Density at top follows Neumann condition
2515	IF ( ocean ) THEN
2516	!$acc parallel present( sums_l )
2517	sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
2518	sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
2519	!$acc end parallel
2520	ENDIF
2521
2522	!
2523	!-- Divergence of vertical flux of resolved scale energy and pressure
2524	!-- fluctuations as well as flux of pressure fluctuation itself (68).
2525	!-- First calculate the products, then the divergence.
2526	!-- Calculation is time consuming. Do it only, if profiles shall be plotted.
2527	IF ( hom(nzb+1,2,55,0) /= 0.0 .OR. hom(nzb+1,2,68,0) /= 0.0 ) THEN
2528
2529	STOP '+++ openACC porting for vertical flux div of resolved scale TKE in flow_statistics is still missing'
2530	sums_ll = 0.0 ! local array
2531
2532	!$OMP DO
2533	DO i = nxl, nxr
2534	DO j = nys, nyn
2535	DO k = nzb_s_inner(j,i)+1, nzt
2536
2537	sums_ll(k,1) = sums_ll(k,1) + 0.5 * w(k,j,i) * ( &
2538	( 0.25 * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1) &
2539	- 0.5 * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) ) &
2540	) )**2 &
2541	+ ( 0.25 * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i) &
2542	- 0.5 * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) ) &
2543	) )**2 &
2544	+ w(k,j,i)**2 )
2545
2546	sums_ll(k,2) = sums_ll(k,2) + 0.5 * w(k,j,i) &
2547	* ( p(k,j,i) + p(k+1,j,i) )
2548
2549	ENDDO
2550	ENDDO
2551	ENDDO
2552	sums_ll(0,1) = 0.0 ! because w is zero at the bottom
2553	sums_ll(nzt+1,1) = 0.0
2554	sums_ll(0,2) = 0.0
2555	sums_ll(nzt+1,2) = 0.0
2556
2557	DO k = nzb+1, nzt
2558	sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
2559	sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
2560	sums_l(k,68,tn) = sums_ll(k,2)
2561	ENDDO
2562	sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
2563	sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
2564	sums_l(nzb,68,tn) = 0.0 ! because w* = 0 at nzb
2565
2566	ENDIF
2567
2568	!
2569	!-- Divergence of vertical flux of SGS TKE and the flux itself (69)
2570	IF ( hom(nzb+1,2,57,0) /= 0.0 .OR. hom(nzb+1,2,69,0) /= 0.0 ) THEN
2571
2572	STOP '+++ openACC porting for vertical flux div of SGS TKE in flow_statistics is still missing'
2573	!$OMP DO
2574	DO i = nxl, nxr
2575	DO j = nys, nyn
2576	DO k = nzb_s_inner(j,i)+1, nzt
2577
2578	sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5 * ( &
2579	(km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
2580	- (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k) &
2581	) * ddzw(k)
2582
2583	sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5 * ( &
2584	(km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
2585	)
2586
2587	ENDDO
2588	ENDDO
2589	ENDDO
2590	sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
2591	sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
2592
2593	ENDIF
2594
2595	!
2596	!-- Horizontal heat fluxes (subgrid, resolved, total).
2597	!-- Do it only, if profiles shall be plotted.
2598	IF ( hom(nzb+1,2,58,0) /= 0.0 ) THEN
2599
2600	STOP '+++ openACC porting for horizontal flux calculation in flow_statistics is still missing'
2601	!$OMP DO
2602	DO i = nxl, nxr
2603	DO j = nys, nyn
2604	DO k = nzb_s_inner(j,i)+1, nzt
2605	!
2606	!-- Subgrid horizontal heat fluxes u"pt", v"pt"
2607	sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5 * &
2608	( kh(k,j,i) + kh(k,j,i-1) ) &
2609	* ( pt(k,j,i-1) - pt(k,j,i) ) &
2610	* ddx * rmask(j,i,sr)
2611	sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5 * &
2612	( kh(k,j,i) + kh(k,j-1,i) ) &
2613	* ( pt(k,j-1,i) - pt(k,j,i) ) &
2614	* ddy * rmask(j,i,sr)
2615	!
2616	!-- Resolved horizontal heat fluxes upt, vpt
2617	sums_l(k,59,tn) = sums_l(k,59,tn) + &
2618	( u(k,j,i) - hom(k,1,1,sr) ) &
2619	* 0.5 * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
2620	pt(k,j,i) - hom(k,1,4,sr) )
2621	pts = 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
2622	pt(k,j,i) - hom(k,1,4,sr) )
2623	sums_l(k,62,tn) = sums_l(k,62,tn) + &
2624	( v(k,j,i) - hom(k,1,2,sr) ) &
2625	* 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
2626	pt(k,j,i) - hom(k,1,4,sr) )
2627	ENDDO
2628	ENDDO
2629	ENDDO
2630	!
2631	!-- Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
2632	sums_l(nzb,58,tn) = 0.0
2633	sums_l(nzb,59,tn) = 0.0
2634	sums_l(nzb,60,tn) = 0.0
2635	sums_l(nzb,61,tn) = 0.0
2636	sums_l(nzb,62,tn) = 0.0
2637	sums_l(nzb,63,tn) = 0.0
2638
2639	ENDIF
2640
2641	!
2642	!-- Calculate the user-defined profiles
2643	CALL user_statistics( 'profiles', sr, tn )
2644	!$OMP END PARALLEL
2645
2646	!
2647	!-- Summation of thread sums
2648	IF ( threads_per_task > 1 ) THEN
2649	STOP '+++ openACC porting for threads_per_task > 1 in flow_statistics is still missing'
2650	DO i = 1, threads_per_task-1
2651	sums_l(:,3,0) = sums_l(:,3,0) + sums_l(:,3,i)
2652	sums_l(:,4:40,0) = sums_l(:,4:40,0) + sums_l(:,4:40,i)
2653	sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
2654	sums_l(:,45:pr_palm,i)
2655	IF ( max_pr_user > 0 ) THEN
2656	sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
2657	sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
2658	sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
2659	ENDIF
2660	ENDDO
2661	ENDIF
2662
2663	!$acc update host( hom, sums, sums_l )
2664
2665	#if defined( __parallel )
2666
2667	!
2668	!-- Compute total sum from local sums
2669	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
2670	CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, &
2671	MPI_SUM, comm2d, ierr )
2672	#else
2673	sums = sums_l(:,:,0)
2674	#endif
2675
2676	!
2677	!-- Final values are obtained by division by the total number of grid points
2678	!-- used for summation. After that store profiles.
2679	!-- Profiles:
2680	DO k = nzb, nzt+1
2681	sums(k,3) = sums(k,3) / ngp_2dh(sr)
2682	sums(k,8:11) = sums(k,8:11) / ngp_2dh_s_inner(k,sr)
2683	sums(k,12:22) = sums(k,12:22) / ngp_2dh(sr)
2684	sums(k,23:29) = sums(k,23:29) / ngp_2dh_s_inner(k,sr)
2685	sums(k,30:32) = sums(k,30:32) / ngp_2dh(sr)
2686	sums(k,33:34) = sums(k,33:34) / ngp_2dh_s_inner(k,sr)
2687	sums(k,35:39) = sums(k,35:39) / ngp_2dh(sr)
2688	sums(k,40) = sums(k,40) / ngp_2dh_s_inner(k,sr)
2689	sums(k,45:53) = sums(k,45:53) / ngp_2dh(sr)
2690	sums(k,54) = sums(k,54) / ngp_2dh_s_inner(k,sr)
2691	sums(k,55:63) = sums(k,55:63) / ngp_2dh(sr)
2692	sums(k,64) = sums(k,64) / ngp_2dh_s_inner(k,sr)
2693	sums(k,65:69) = sums(k,65:69) / ngp_2dh(sr)
2694	sums(k,70:pr_palm-2) = sums(k,70:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
2695	ENDDO
2696
2697	!-- Upstream-parts
2698	sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
2699	!-- u* and so on
2700	!-- As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
2701	!-- size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
2702	!-- above the topography, they are being divided by ngp_2dh(sr)
2703	sums(nzb:nzb+3,pr_palm) = sums(nzb:nzb+3,pr_palm) / &
2704	ngp_2dh(sr)
2705	sums(nzb+12,pr_palm) = sums(nzb+12,pr_palm) / & ! qs
2706	ngp_2dh(sr)
2707	!-- eges, e*
2708	sums(nzb+4:nzb+5,pr_palm) = sums(nzb+4:nzb+5,pr_palm) / &
2709	ngp_3d(sr)
2710	!-- Old and new divergence
2711	sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
2712	ngp_3d_inner(sr)
2713
2714	!-- User-defined profiles
2715	IF ( max_pr_user > 0 ) THEN
2716	DO k = nzb, nzt+1
2717	sums(k,pr_palm+1:pr_palm+max_pr_user) = &
2718	sums(k,pr_palm+1:pr_palm+max_pr_user) / &
2719	ngp_2dh_s_inner(k,sr)
2720	ENDDO
2721	ENDIF
2722
2723	!
2724	!-- Collect horizontal average in hom.
2725	!-- Compute deduced averages (e.g. total heat flux)
2726	hom(:,1,3,sr) = sums(:,3) ! w
2727	hom(:,1,8,sr) = sums(:,8) ! e profiles 5-7 are initial profiles
2728	hom(:,1,9,sr) = sums(:,9) ! km
2729	hom(:,1,10,sr) = sums(:,10) ! kh
2730	hom(:,1,11,sr) = sums(:,11) ! l
2731	hom(:,1,12,sr) = sums(:,12) ! w"u"
2732	hom(:,1,13,sr) = sums(:,13) ! wu
2733	hom(:,1,14,sr) = sums(:,14) ! w"v"
2734	hom(:,1,15,sr) = sums(:,15) ! wv
2735	hom(:,1,16,sr) = sums(:,16) ! w"pt"
2736	hom(:,1,17,sr) = sums(:,17) ! wpt
2737	hom(:,1,18,sr) = sums(:,16) + sums(:,17) ! wpt
2738	hom(:,1,19,sr) = sums(:,12) + sums(:,13) ! wu
2739	hom(:,1,20,sr) = sums(:,14) + sums(:,15) ! wv
2740	hom(:,1,21,sr) = sums(:,21) ! wptBC
2741	hom(:,1,22,sr) = sums(:,16) + sums(:,21) ! wptBC
2742	! profile 24 is initial profile (sa)
2743	! profiles 25-29 left empty for initial
2744	! profiles
2745	hom(:,1,30,sr) = sums(:,30) ! u*2
2746	hom(:,1,31,sr) = sums(:,31) ! v*2
2747	hom(:,1,32,sr) = sums(:,32) ! w*2
2748	hom(:,1,33,sr) = sums(:,33) ! pt*2
2749	hom(:,1,34,sr) = sums(:,34) ! e*
2750	hom(:,1,35,sr) = sums(:,35) ! w2pt
2751	hom(:,1,36,sr) = sums(:,36) ! wpt2
2752	hom(:,1,37,sr) = sums(:,37) ! we
2753	hom(:,1,38,sr) = sums(:,38) ! w*3
2754	hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20 )**1.5 ! Sw
2755	hom(:,1,40,sr) = sums(:,40) ! p
2756	hom(:,1,45,sr) = sums(:,45) ! w"vpt"
2757	hom(:,1,46,sr) = sums(:,46) ! wvpt
2758	hom(:,1,47,sr) = sums(:,45) + sums(:,46) ! wvpt
2759	hom(:,1,48,sr) = sums(:,48) ! w"q" (w"qv")
2760	hom(:,1,49,sr) = sums(:,49) ! wq (wqv)
2761	hom(:,1,50,sr) = sums(:,48) + sums(:,49) ! wq (wqv)
2762	hom(:,1,51,sr) = sums(:,51) ! w"qv"
2763	hom(:,1,52,sr) = sums(:,52) ! wqv
2764	hom(:,1,53,sr) = sums(:,52) + sums(:,51) ! wq (wqv)
2765	hom(:,1,54,sr) = sums(:,54) ! ql
2766	hom(:,1,55,sr) = sums(:,55) ! wuu*/dz
2767	hom(:,1,56,sr) = sums(:,56) ! wp/dz
2768	hom(:,1,57,sr) = sums(:,57) ! ( w"e + w"p"/rho )/dz
2769	hom(:,1,58,sr) = sums(:,58) ! u"pt"
2770	hom(:,1,59,sr) = sums(:,59) ! upt
2771	hom(:,1,60,sr) = sums(:,58) + sums(:,59) ! upt_t
2772	hom(:,1,61,sr) = sums(:,61) ! v"pt"
2773	hom(:,1,62,sr) = sums(:,62) ! vpt
2774	hom(:,1,63,sr) = sums(:,61) + sums(:,62) ! vpt_t
2775	hom(:,1,64,sr) = sums(:,64) ! rho
2776	hom(:,1,65,sr) = sums(:,65) ! w"sa"
2777	hom(:,1,66,sr) = sums(:,66) ! wsa
2778	hom(:,1,67,sr) = sums(:,65) + sums(:,66) ! wsa
2779	hom(:,1,68,sr) = sums(:,68) ! wp
2780	hom(:,1,69,sr) = sums(:,69) ! w"e + w"p"/rho
2781	hom(:,1,70,sr) = sums(:,70) ! q*2
2782	hom(:,1,71,sr) = sums(:,71) ! prho
2783	hom(:,1,72,sr) = hyp * 1E-4 ! hyp in dbar
2784	hom(:,1,73,sr) = sums(:,73) ! nr
2785	hom(:,1,74,sr) = sums(:,74) ! qr
2786	hom(:,1,75,sr) = sums(:,75) ! qc
2787	hom(:,1,76,sr) = sums(:,76) ! prr (precipitation rate)
2788	! 77 is initial density profile
2789	hom(:,1,78,sr) = ug ! ug
2790	hom(:,1,79,sr) = vg ! vg
2791	hom(:,1,80,sr) = w_subs ! w_subs
2792
2793	hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
2794	! upstream-parts u_x, u_y, u_z, v_x,
2795	! v_y, usw. (in last but one profile)
2796	hom(:,1,pr_palm,sr) = sums(:,pr_palm)
2797	! u, w'u', w'v', t (in last profile)
2798
2799	IF ( max_pr_user > 0 ) THEN ! user-defined profiles
2800	hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
2801	sums(:,pr_palm+1:pr_palm+max_pr_user)
2802	ENDIF
2803
2804	!
2805	!-- Determine the boundary layer height using two different schemes.
2806	!-- First scheme: Starting from the Earth's (Ocean's) surface, look for the
2807	!-- first relative minimum (maximum) of the total heat flux.
2808	!-- The corresponding height is assumed as the boundary layer height, if it
2809	!-- is less than 1.5 times the height where the heat flux becomes negative
2810	!-- (positive) for the first time.
2811	z_i(1) = 0.0
2812	first = .TRUE.
2813
2814	IF ( ocean ) THEN
2815	DO k = nzt, nzb+1, -1
2816	IF ( first .AND. hom(k,1,18,sr) < 0.0 &
2817	.AND. abs(hom(k,1,18,sr)) > 1.0E-8) THEN
2818	first = .FALSE.
2819	height = zw(k)
2820	ENDIF
2821	IF ( hom(k,1,18,sr) < 0.0 .AND. &
2822	abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
2823	hom(k-1,1,18,sr) > hom(k,1,18,sr) ) THEN
2824	IF ( zw(k) < 1.5 * height ) THEN
2825	z_i(1) = zw(k)
2826	ELSE
2827	z_i(1) = height
2828	ENDIF
2829	EXIT
2830	ENDIF
2831	ENDDO
2832	ELSE
2833	DO k = nzb, nzt-1
2834	IF ( first .AND. hom(k,1,18,sr) < 0.0 &
2835	.AND. abs(hom(k,1,18,sr)) > 1.0E-8 ) THEN
2836	first = .FALSE.
2837	height = zw(k)
2838	ENDIF
2839	IF ( hom(k,1,18,sr) < 0.0 .AND. &
2840	abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
2841	hom(k+1,1,18,sr) > hom(k,1,18,sr) ) THEN
2842	IF ( zw(k) < 1.5 * height ) THEN
2843	z_i(1) = zw(k)
2844	ELSE
2845	z_i(1) = height
2846	ENDIF
2847	EXIT
2848	ENDIF
2849	ENDDO
2850	ENDIF
2851
2852	!
2853	!-- Second scheme: Gradient scheme from Sullivan et al. (1998), modified
2854	!-- by Uhlenbrock(2006). The boundary layer height is the height with the
2855	!-- maximal local temperature gradient: starting from the second (the last
2856	!-- but one) vertical gridpoint, the local gradient must be at least
2857	!-- 0.2K/100m and greater than the next four gradients.
2858	!-- WARNING: The threshold value of 0.2K/100m must be adjusted for the
2859	!-- ocean case!
2860	z_i(2) = 0.0
2861	DO k = nzb+1, nzt+1
2862	dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
2863	ENDDO
2864	dptdz_threshold = 0.2_wp / 100.0_wp
2865
2866	IF ( ocean ) THEN
2867	DO k = nzt+1, nzb+5, -1
2868	IF ( dptdz(k) > dptdz_threshold .AND. &
2869	dptdz(k) > dptdz(k-1) .AND. dptdz(k) > dptdz(k-2) .AND. &
2870	dptdz(k) > dptdz(k-3) .AND. dptdz(k) > dptdz(k-4) ) THEN
2871	z_i(2) = zw(k-1)
2872	EXIT
2873	ENDIF
2874	ENDDO
2875	ELSE
2876	DO k = nzb+1, nzt-3
2877	IF ( dptdz(k) > dptdz_threshold .AND. &
2878	dptdz(k) > dptdz(k+1) .AND. dptdz(k) > dptdz(k+2) .AND. &
2879	dptdz(k) > dptdz(k+3) .AND. dptdz(k) > dptdz(k+4) ) THEN
2880	z_i(2) = zw(k-1)
2881	EXIT
2882	ENDIF
2883	ENDDO
2884	ENDIF
2885
2886	hom(nzb+6,1,pr_palm,sr) = z_i(1)
2887	hom(nzb+7,1,pr_palm,sr) = z_i(2)
2888
2889	!
2890	!-- Computation of both the characteristic vertical velocity and
2891	!-- the characteristic convective boundary layer temperature.
2892	!-- The horizontal average at nzb+1 is input for the average temperature.
2893	IF ( hom(nzb,1,18,sr) > 0.0 .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8 &
2894	.AND. z_i(1) /= 0.0 ) THEN
2895	hom(nzb+8,1,pr_palm,sr) = ( g / hom(nzb+1,1,4,sr) * &
2896	hom(nzb,1,18,sr) * &
2897	ABS( z_i(1) ) )**0.333333333
2898	!-- so far this only works if Prandtl layer is used
2899	hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
2900	ELSE
2901	hom(nzb+8,1,pr_palm,sr) = 0.0
2902	hom(nzb+11,1,pr_palm,sr) = 0.0
2903	ENDIF
2904
2905	!
2906	!-- Collect the time series quantities
2907	ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr) ! E
2908	ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr) ! E*
2909	ts_value(3,sr) = dt_3d
2910	ts_value(4,sr) = hom(nzb,1,pr_palm,sr) ! u*
2911	ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr) ! th*
2912	ts_value(6,sr) = u_max
2913	ts_value(7,sr) = v_max
2914	ts_value(8,sr) = w_max
2915	ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr) ! new divergence
2916	ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr) ! old Divergence
2917	ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr) ! z_i(1)
2918	ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr) ! z_i(2)
2919	ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr) ! w*
2920	ts_value(14,sr) = hom(nzb,1,16,sr) ! w'pt' at k=0
2921	ts_value(15,sr) = hom(nzb+1,1,16,sr) ! w'pt' at k=1
2922	ts_value(16,sr) = hom(nzb+1,1,18,sr) ! wpt at k=1
2923	ts_value(17,sr) = hom(nzb,1,4,sr) ! pt(0)
2924	ts_value(18,sr) = hom(nzb+1,1,4,sr) ! pt(zp)
2925	ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr) ! u'w' at k=0
2926	ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr) ! v'w' at k=0
2927	ts_value(21,sr) = hom(nzb,1,48,sr) ! w"q" at k=0
2928
2929	IF ( ts_value(5,sr) /= 0.0 ) THEN
2930	ts_value(22,sr) = ts_value(4,sr)**2 / &
2931	( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
2932	ELSE
2933	ts_value(22,sr) = 10000.0
2934	ENDIF
2935
2936	ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr) ! q*
2937	!
2938	!-- Calculate additional statistics provided by the user interface
2939	CALL user_statistics( 'time_series', sr, 0 )
2940
2941	ENDDO ! loop of the subregions
2942
2943	!$acc end data
2944
2945	!
2946	!-- If required, sum up horizontal averages for subsequent time averaging
2947	IF ( do_sum ) THEN
2948	IF ( average_count_pr == 0 ) hom_sum = 0.0
2949	hom_sum = hom_sum + hom(:,1,:,:)
2950	average_count_pr = average_count_pr + 1
2951	do_sum = .FALSE.
2952	ENDIF
2953
2954	!
2955	!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
2956	!-- This flag is reset after each time step in time_integration.
2957	flow_statistics_called = .TRUE.
2958
2959	CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
2960
2961
2962	END SUBROUTINE flow_statistics
2963	#endif

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |