Home

Context Navigation

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

Last change on this file since 1817 was 1816, checked in by raasch, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 147.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |