Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |