Home

Context Navigation

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

Last change on this file since 1682 was 1682, checked in by knoop, 9 years ago
Code annotations made doxygen readable
Property svn:keywords set to `Id`
File size: 144.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |