Home

Context Navigation

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

Last change on this file since 1396 was 1396, checked in by raasch, 10 years ago
bugfix: "copyin" replaced by "update device" in openacc-branch
Property svn:keywords set to `Id`
File size: 131.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |