Home

Context Navigation

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

Last change on this file since 1567 was 1567, checked in by suehring, 9 years ago
Bugfixes in monotonic limter.
Property svn:keywords set to `Id`
File size: 138.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |