Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |