Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |