Home

Context Navigation

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

Last change on this file since 1498 was 1498, checked in by suehring, 9 years ago
Comments added
Property svn:keywords set to `Id`
File size: 133.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |