Home

Context Navigation

flow_statistics.f90 @ 1257

Last change on this file since 1257 was 1257, checked in by raasch, 10 years ago

New:
---

openACC porting of timestep calculation
(modules, timestep, time_integration)

Changed:

openACC loop directives and vector clauses removed (because they do not give any performance improvement with PGI
compiler versions > 13.6)
(advec_ws, buoyancy, coriolis, diffusion_e, diffusion_s, diffusion_u, diffusion_v, diffusion_w, diffusivities, exchange_horiz, fft_xy, pres, production_e, transpose, tridia_solver, wall_fluxes)

openACC loop independent clauses added
(boundary_conds, prandtl_fluxes, pres)

openACC declare create statements moved after FORTRAN declaration statement
(diffusion_u, diffusion_v, diffusion_w, fft_xy, poisfft, production_e, tridia_solver)

openACC end parallel replaced by end parallel loop
(flow_statistics, pres)

openACC "kernels do" replaced by "kernels loop"
(prandtl_fluxes)

output format for theta* changed to avoid output of *
(run_control)

Errors:

bugfix for calculation of advective timestep (old version may cause wrong timesteps in case of
vertixcally stretched grids)
Attention: standard run-control output has changed!
(timestep)

Property svn:keywords set to Id

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

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

Context Navigation

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

Download in other formats: