Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |