Home

Context Navigation

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

Last change on this file since 1319 was 1319, checked in by raasch, 10 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 120.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |