Home

Context Navigation

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

Last change on this file since 1965 was 1961, checked in by suehring, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 150.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |