Home

Context Navigation

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

Last change on this file since 1853 was 1853, checked in by maronga, 8 years ago
bugfix: radiation_scheme = constant caused model crash when used in combination with land surface model
Property svn:keywords set to `Id`
File size: 147.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |