Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |