Home

Context Navigation

flow_statistics.f90 @ 1221

Last change on this file since 1221 was 1221, checked in by raasch, 11 years ago

New:

openACC porting of reduction operations
additional 3D-flag arrays for replacing the 2D-index arrays nzb_s_inner and nzb_diff_s_inner
(flow_statistics, init_grid, init_3d_model, modules, palm, pres, time_integration)

Changed:

for PGI/openACC performance reasons set default compile options for openACC to "-ta=nocache",
and set environment variable PGI_ACC_SYNCHRONOUS=1
(MAKE.inc.pgi.openacc, palm_simple_run)

wall_flags_0 changed to 32bit INTEGER, additional array wall_flags_00 introduced to hold
bits 32-63
(advec_ws, init_grid, modules, palm)

Errors:

dummy argument tri in 1d-routines replaced by tri_for_1d because of name
conflict with arry tri in module arrays_3d
(tridia_solver)

Property svn:keywords set to Id

File size: 120.0 KB

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

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

Context Navigation

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

Download in other formats: