Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |