Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |