Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |