Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |