Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |