Home

Context Navigation

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

Last change on this file since 2030 was 2027, checked in by suehring, 8 years ago
Last commit documented
Property svn:keywords set to `Id`
File size: 150.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |