Home

Context Navigation

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

Last change on this file since 2000 was 2000, checked in by knoop, 8 years ago
Forced header and separation lines into 80 columns
Property svn:keywords set to `Id`
File size: 150.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |