Home

Context Navigation

flow_statistics.f90 @ 1179

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

New:
---
Initial profiles can be used as reference state in the buoyancy term. New parameter
reference_state introduced. Calculation and handling of reference state in buoyancy term revised.
binary version for restart files changed from 3.9 to 3.9a (no downward compatibility!),
initial profile for rho added to hom (id=77)

Errors:

small bugfix for background communication (time_integration)

Property svn:keywords set to Id

File size: 56.6 KB

Rev	Line
[1]	1	SUBROUTINE flow_statistics
	2
[1036]	3	!--------------------------------------------------------------------------------!
	4	! This file is part of PALM.
	5	!
	6	! PALM is free software: you can redistribute it and/or modify it under the terms
	7	! of the GNU General Public License as published by the Free Software Foundation,
	8	! either version 3 of the License, or (at your option) any later version.
	9	!
	10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	13	!
	14	! You should have received a copy of the GNU General Public License along with
	15	! PALM. If not, see <http://www.gnu.org/licenses/>.
	16	!
	17	! Copyright 1997-2012 Leibniz University Hannover
	18	!--------------------------------------------------------------------------------!
	19	!
[254]	20	! Current revisions:
[1]	21	! -----------------
[1179]	22	! comment for profile 77 added
[1008]	23	!
	24	! Former revisions:
	25	! -----------------
	26	! $Id: flow_statistics.f90 1179 2013-06-14 05:57:58Z raasch $
	27	!
[1116]	28	! 1115 2013-03-26 18:16:16Z hoffmann
	29	! ql is calculated by calc_liquid_water_content
	30	!
[1112]	31	! 1111 2013-03-08 23:54:10Z raasch
	32	! openACC directive added
	33	!
[1054]	34	! 1053 2012-11-13 17:11:03Z hoffmann
[1112]	35	! additions for two-moment cloud physics scheme:
[1054]	36	! +nr, qr, qc, prr
	37	!
[1037]	38	! 1036 2012-10-22 13:43:42Z raasch
	39	! code put under GPL (PALM 3.9)
	40	!
[1008]	41	! 1007 2012-09-19 14:30:36Z franke
[1007]	42	! Calculation of buoyancy flux for humidity in case of WS-scheme is now using
	43	! turbulent fluxes of WS-scheme
	44	! Bugfix: Calculation of subgridscale buoyancy flux for humidity and cloud
	45	! droplets at nzb and nzt added
[700]	46	!
[802]	47	! 801 2012-01-10 17:30:36Z suehring
	48	! Calculation of turbulent fluxes in advec_ws is now thread-safe.
	49	!
[744]	50	! 743 2011-08-18 16:10:16Z suehring
	51	! Calculation of turbulent fluxes with WS-scheme only for the whole model
	52	! domain, not for user-defined subregions.
	53	!
[710]	54	! 709 2011-03-30 09:31:40Z raasch
	55	! formatting adjustments
	56	!
[700]	57	! 699 2011-03-22 17:52:22Z suehring
[699]	58	! Bugfix in calculation of vertical velocity skewness. The added absolute value
	59	! avoid negative values in the root. Negative values of w'w' can occur at the
	60	! top or bottom of the model domain due to degrading the order of advection
	61	! scheme. Furthermore the calculation will be the same for all advection
	62	! schemes.
[392]	63	!
[697]	64	! 696 2011-03-18 07:03:49Z raasch
	65	! Bugfix: Summation of Wicker-Skamarock scheme fluxes and variances for all
	66	! threads
	67	!
[679]	68	! 678 2011-02-02 14:31:56Z raasch
	69	! Bugfix in calculation of the divergence of vertical flux of resolved scale
	70	! energy, pressure fluctuations, and flux of pressure fluctuation itself
	71	!
[674]	72	! 673 2011-01-18 16:19:48Z suehring
	73	! Top bc for the horizontal velocity variances added for ocean runs.
	74	! Setting the corresponding bottom bc moved to advec_ws.
	75	!
[668]	76	! 667 2010-12-23 12:06:00Z suehring/gryschka
	77	! When advection is computed with ws-scheme, turbulent fluxes are already
	78	! computed in the respective advection routines and buffered in arrays
	79	! sums_xx_ws_l(). This is due to a consistent treatment of statistics with the
	80	! numerics and to avoid unphysical kinks near the surface.
	81	! So some if requests has to be done to dicern between fluxes from ws-scheme
	82	! other advection schemes.
	83	! Furthermore the computation of z_i is only done if the heat flux exceeds a
	84	! minimum value. This affects only simulations of a neutral boundary layer and
	85	! is due to reasons of computations in the advection scheme.
	86	!
[625]	87	! 624 2010-12-10 11:46:30Z heinze
	88	! Calculation of q*2 added
	89	!
[623]	90	! 622 2010-12-10 08:08:13Z raasch
	91	! optional barriers included in order to speed up collective operations
	92	!
[392]	93	! 388 2009-09-23 09:40:33Z raasch
[388]	94	! Vertical profiles of potential density and hydrostatic pressure are
	95	! calculated.
[343]	96	! Added missing timeseries calculation of w"q"(0), moved timeseries q* to the
	97	! end.
[291]	98	! Temperature gradient criterion for estimating the boundary layer height
	99	! replaced by the gradient criterion of Sullivan et al. (1998).
[254]	100	! Output of messages replaced by message handling routine.
[1]	101	!
[198]	102	! 197 2008-09-16 15:29:03Z raasch
	103	! Spline timeseries splptx etc. removed, timeseries w'u', w'v', w'q' (k=0)
	104	! added,
	105	! bugfix: divide sums(k,8) (e) and sums(k,34) (e*) by ngp_2dh_s_inner(k,sr)
	106	! (like other scalars)
	107	!
[139]	108	! 133 2007-11-20 10:10:53Z letzel
	109	! Vertical profiles now based on nzb_s_inner; they are divided by
	110	! ngp_2dh_s_inner (scalars, procucts of scalars) and ngp_2dh (staggered
	111	! velocity components and their products, procucts of scalars and velocity
	112	! components), respectively.
	113	!
[110]	114	! 106 2007-08-16 14:30:26Z raasch
	115	! Prescribed momentum fluxes at the top surface are used,
	116	! profiles for wp and w"e are calculated
	117	!
[98]	118	! 97 2007-06-21 08:23:15Z raasch
	119	! Statistics for ocean version (salinity, density) added,
	120	! calculation of z_i and Deardorff velocity scale adjusted to be used with
	121	! the ocean version
	122	!
[90]	123	! 87 2007-05-22 15:46:47Z raasch
	124	! Two more arguments added to user_statistics, which is now also called for
	125	! user-defined profiles,
	126	! var_hom and var_sum renamed pr_palm
	127	!
[83]	128	! 82 2007-04-16 15:40:52Z raasch
	129	! Cpp-directive lcmuk changed to intel_openmp_bug
	130	!
[77]	131	! 75 2007-03-22 09:54:05Z raasch
	132	! Collection of time series quantities moved from routine flow_statistics to
	133	! here, routine user_statistics is called for each statistic region,
	134	! moisture renamed humidity
	135	!
[39]	136	! 19 2007-02-23 04:53:48Z raasch
[77]	137	! fluxes at top modified (tswst, qswst)
[39]	138	!
[3]	139	! RCS Log replace by Id keyword, revision history cleaned up
	140	!
[1]	141	! Revision 1.41 2006/08/04 14:37:50 raasch
	142	! Error removed in non-parallel part (sums_l)
	143	!
	144	! Revision 1.1 1997/08/11 06:15:17 raasch
	145	! Initial revision
	146	!
	147	!
	148	! Description:
	149	! ------------
	150	! Compute average profiles and further average flow quantities for the different
	151	! user-defined (sub-)regions. The region indexed 0 is the total model domain.
	152	!
[132]	153	! NOTE: For simplicity, nzb_s_inner and nzb_diff_s_inner are being used as a
	154	! ---- lower vertical index for k-loops for all variables, although strictly
	155	! speaking the k-loops would have to be split up according to the staggered
	156	! grid. However, this implies no error since staggered velocity components are
	157	! zero at the walls and inside buildings.
[1]	158	!------------------------------------------------------------------------------!
	159
	160	USE arrays_3d
	161	USE cloud_parameters
[709]	162	USE control_parameters
[1]	163	USE cpulog
	164	USE grid_variables
	165	USE indices
	166	USE interfaces
	167	USE pegrid
	168	USE statistics
	169
	170	IMPLICIT NONE
	171
	172	INTEGER :: i, j, k, omp_get_thread_num, sr, tn
	173	LOGICAL :: first
[291]	174	REAL :: dptdz_threshold, height, pts, sums_l_eper, sums_l_etot, ust, &
	175	ust2, u2, vst, vst2, v2, w2, z_i(2)
	176	REAL :: dptdz(nzb+1:nzt+1)
[1]	177	REAL :: sums_ll(nzb:nzt+1,2)
	178
	179	CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
	180
	181	!
	182	!-- To be on the safe side, check whether flow_statistics has already been
	183	!-- called once after the current time step
	184	IF ( flow_statistics_called ) THEN
[254]	185
[274]	186	message_string = 'flow_statistics is called two times within one ' // &
	187	'timestep'
[254]	188	CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
[1007]	189
[1]	190	ENDIF
	191
[1111]	192	!$acc update host( km, kh, e, pt, qs, qsws, rif, shf, ts, u, v, w )
	193
[1]	194	!
	195	!-- Compute statistics for each (sub-)region
	196	DO sr = 0, statistic_regions
	197
	198	!
	199	!-- Initialize (local) summation array
	200	sums_l = 0.0
	201
	202	!
	203	!-- Store sums that have been computed in other subroutines in summation
	204	!-- array
	205	sums_l(:,11,:) = sums_l_l(:,sr,:) ! mixing length from diffusivities
	206	!-- WARNING: next line still has to be adjusted for OpenMP
	207	sums_l(:,21,0) = sums_wsts_bc_l(:,sr) ! heat flux from advec_s_bc
[87]	208	sums_l(nzb+9,pr_palm,0) = sums_divold_l(sr) ! old divergence from pres
	209	sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr) ! new divergence from pres
[1]	210
[667]	211	!
	212	!-- Copy the turbulent quantities, evaluated in the advection routines to
	213	!-- the local array sums_l() for further computations
[743]	214	IF ( ws_scheme_mom .AND. sr == 0 ) THEN
[696]	215
[1007]	216	!
[673]	217	!-- According to the Neumann bc for the horizontal velocity components,
	218	!-- the corresponding fluxes has to satisfiy the same bc.
	219	IF ( ocean ) THEN
[801]	220	sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
[1007]	221	sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)
[673]	222	ENDIF
[696]	223
	224	DO i = 0, threads_per_task-1
[1007]	225	!
[696]	226	!-- Swap the turbulent quantities evaluated in advec_ws.
[801]	227	sums_l(:,13,i) = sums_wsus_ws_l(:,i) ! wu
	228	sums_l(:,15,i) = sums_wsvs_ws_l(:,i) ! wv
	229	sums_l(:,30,i) = sums_us2_ws_l(:,i) ! u*2
	230	sums_l(:,31,i) = sums_vs2_ws_l(:,i) ! v*2
	231	sums_l(:,32,i) = sums_ws2_ws_l(:,i) ! w*2
[696]	232	sums_l(:,34,i) = sums_l(:,34,i) + 0.5 * &
[801]	233	( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) + &
	234	sums_ws2_ws_l(:,i) ) ! e*
[696]	235	DO k = nzb, nzt
[801]	236	sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5 * ( &
	237	sums_us2_ws_l(k,i) + &
	238	sums_vs2_ws_l(k,i) + &
	239	sums_ws2_ws_l(k,i) )
[696]	240	ENDDO
[667]	241	ENDDO
[696]	242
[667]	243	ENDIF
[696]	244
[743]	245	IF ( ws_scheme_sca .AND. sr == 0 ) THEN
[696]	246
	247	DO i = 0, threads_per_task-1
[801]	248	sums_l(:,17,i) = sums_wspts_ws_l(:,i) ! wpt from advec_s_ws
	249	IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! wsa
[696]	250	IF ( humidity .OR. passive_scalar ) sums_l(:,49,i) = &
[801]	251	sums_wsqs_ws_l(:,i) !wq
[696]	252	ENDDO
	253
[667]	254	ENDIF
[305]	255	!
[1]	256	!-- Horizontally averaged profiles of horizontal velocities and temperature.
	257	!-- They must have been computed before, because they are already required
	258	!-- for other horizontal averages.
	259	tn = 0
[667]	260
[1]	261	!$OMP PARALLEL PRIVATE( i, j, k, tn )
[82]	262	#if defined( __intel_openmp_bug )
[1]	263	tn = omp_get_thread_num()
	264	#else
	265	!$ tn = omp_get_thread_num()
	266	#endif
	267
	268	!$OMP DO
	269	DO i = nxl, nxr
	270	DO j = nys, nyn
[132]	271	DO k = nzb_s_inner(j,i), nzt+1
[1]	272	sums_l(k,1,tn) = sums_l(k,1,tn) + u(k,j,i) * rmask(j,i,sr)
	273	sums_l(k,2,tn) = sums_l(k,2,tn) + v(k,j,i) * rmask(j,i,sr)
	274	sums_l(k,4,tn) = sums_l(k,4,tn) + pt(k,j,i) * rmask(j,i,sr)
	275	ENDDO
	276	ENDDO
	277	ENDDO
	278
	279	!
[96]	280	!-- Horizontally averaged profile of salinity
	281	IF ( ocean ) THEN
	282	!$OMP DO
	283	DO i = nxl, nxr
	284	DO j = nys, nyn
[132]	285	DO k = nzb_s_inner(j,i), nzt+1
[96]	286	sums_l(k,23,tn) = sums_l(k,23,tn) + &
	287	sa(k,j,i) * rmask(j,i,sr)
	288	ENDDO
	289	ENDDO
	290	ENDDO
	291	ENDIF
	292
	293	!
[1]	294	!-- Horizontally averaged profiles of virtual potential temperature,
	295	!-- total water content, specific humidity and liquid water potential
	296	!-- temperature
[75]	297	IF ( humidity ) THEN
[1]	298	!$OMP DO
	299	DO i = nxl, nxr
	300	DO j = nys, nyn
[132]	301	DO k = nzb_s_inner(j,i), nzt+1
[1]	302	sums_l(k,44,tn) = sums_l(k,44,tn) + &
	303	vpt(k,j,i) * rmask(j,i,sr)
	304	sums_l(k,41,tn) = sums_l(k,41,tn) + &
	305	q(k,j,i) * rmask(j,i,sr)
	306	ENDDO
	307	ENDDO
	308	ENDDO
	309	IF ( cloud_physics ) THEN
	310	!$OMP DO
	311	DO i = nxl, nxr
	312	DO j = nys, nyn
[132]	313	DO k = nzb_s_inner(j,i), nzt+1
[1]	314	sums_l(k,42,tn) = sums_l(k,42,tn) + &
	315	( q(k,j,i) - ql(k,j,i) ) * rmask(j,i,sr)
	316	sums_l(k,43,tn) = sums_l(k,43,tn) + ( &
	317	pt(k,j,i) + l_d_cppt_d_t(k) ql(k,j,i) &
	318	) * rmask(j,i,sr)
	319	ENDDO
	320	ENDDO
	321	ENDDO
	322	ENDIF
	323	ENDIF
	324
	325	!
	326	!-- Horizontally averaged profiles of passive scalar
	327	IF ( passive_scalar ) THEN
	328	!$OMP DO
	329	DO i = nxl, nxr
	330	DO j = nys, nyn
[132]	331	DO k = nzb_s_inner(j,i), nzt+1
[1]	332	sums_l(k,41,tn) = sums_l(k,41,tn) + q(k,j,i) * rmask(j,i,sr)
	333	ENDDO
	334	ENDDO
	335	ENDDO
	336	ENDIF
	337	!$OMP END PARALLEL
	338	!
	339	!-- Summation of thread sums
	340	IF ( threads_per_task > 1 ) THEN
	341	DO i = 1, threads_per_task-1
	342	sums_l(:,1,0) = sums_l(:,1,0) + sums_l(:,1,i)
	343	sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i)
	344	sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i)
[96]	345	IF ( ocean ) THEN
	346	sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)
	347	ENDIF
[75]	348	IF ( humidity ) THEN
[1]	349	sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
	350	sums_l(:,44,0) = sums_l(:,44,0) + sums_l(:,44,i)
	351	IF ( cloud_physics ) THEN
	352	sums_l(:,42,0) = sums_l(:,42,0) + sums_l(:,42,i)
	353	sums_l(:,43,0) = sums_l(:,43,0) + sums_l(:,43,i)
	354	ENDIF
	355	ENDIF
	356	IF ( passive_scalar ) THEN
	357	sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
	358	ENDIF
	359	ENDDO
	360	ENDIF
	361
	362	#if defined( __parallel )
	363	!
	364	!-- Compute total sum from local sums
[622]	365	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[1]	366	CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL, &
	367	MPI_SUM, comm2d, ierr )
[622]	368	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[1]	369	CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL, &
	370	MPI_SUM, comm2d, ierr )
[622]	371	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[1]	372	CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL, &
	373	MPI_SUM, comm2d, ierr )
[96]	374	IF ( ocean ) THEN
[622]	375	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[96]	376	CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb, &
	377	MPI_REAL, MPI_SUM, comm2d, ierr )
	378	ENDIF
[75]	379	IF ( humidity ) THEN
[622]	380	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[1]	381	CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb, &
	382	MPI_REAL, MPI_SUM, comm2d, ierr )
[622]	383	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[1]	384	CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
	385	MPI_REAL, MPI_SUM, comm2d, ierr )
	386	IF ( cloud_physics ) THEN
[622]	387	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[1]	388	CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb, &
	389	MPI_REAL, MPI_SUM, comm2d, ierr )
[622]	390	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[1]	391	CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb, &
	392	MPI_REAL, MPI_SUM, comm2d, ierr )
	393	ENDIF
	394	ENDIF
	395
	396	IF ( passive_scalar ) THEN
[622]	397	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[1]	398	CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
	399	MPI_REAL, MPI_SUM, comm2d, ierr )
	400	ENDIF
	401	#else
	402	sums(:,1) = sums_l(:,1,0)
	403	sums(:,2) = sums_l(:,2,0)
	404	sums(:,4) = sums_l(:,4,0)
[96]	405	IF ( ocean ) sums(:,23) = sums_l(:,23,0)
[75]	406	IF ( humidity ) THEN
[1]	407	sums(:,44) = sums_l(:,44,0)
	408	sums(:,41) = sums_l(:,41,0)
	409	IF ( cloud_physics ) THEN
	410	sums(:,42) = sums_l(:,42,0)
	411	sums(:,43) = sums_l(:,43,0)
	412	ENDIF
	413	ENDIF
	414	IF ( passive_scalar ) sums(:,41) = sums_l(:,41,0)
	415	#endif
	416
	417	!
	418	!-- Final values are obtained by division by the total number of grid points
	419	!-- used for summation. After that store profiles.
[132]	420	sums(:,1) = sums(:,1) / ngp_2dh(sr)
	421	sums(:,2) = sums(:,2) / ngp_2dh(sr)
	422	sums(:,4) = sums(:,4) / ngp_2dh_s_inner(:,sr)
[1]	423	hom(:,1,1,sr) = sums(:,1) ! u
	424	hom(:,1,2,sr) = sums(:,2) ! v
	425	hom(:,1,4,sr) = sums(:,4) ! pt
	426
[667]	427
[1]	428	!
[96]	429	!-- Salinity
	430	IF ( ocean ) THEN
[132]	431	sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr)
[96]	432	hom(:,1,23,sr) = sums(:,23) ! sa
	433	ENDIF
	434
	435	!
[1]	436	!-- Humidity and cloud parameters
[75]	437	IF ( humidity ) THEN
[132]	438	sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr)
	439	sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
[1]	440	hom(:,1,44,sr) = sums(:,44) ! vpt
	441	hom(:,1,41,sr) = sums(:,41) ! qv (q)
	442	IF ( cloud_physics ) THEN
[132]	443	sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr)
	444	sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr)
[1]	445	hom(:,1,42,sr) = sums(:,42) ! qv
	446	hom(:,1,43,sr) = sums(:,43) ! pt
	447	ENDIF
	448	ENDIF
	449
	450	!
	451	!-- Passive scalar
[132]	452	IF ( passive_scalar ) hom(:,1,41,sr) = sums(:,41) / &
	453	ngp_2dh_s_inner(:,sr) ! s (q)
[1]	454
	455	!
	456	!-- Horizontally averaged profiles of the remaining prognostic variables,
	457	!-- variances, the total and the perturbation energy (single values in last
	458	!-- column of sums_l) and some diagnostic quantities.
[132]	459	!-- NOTE: for simplicity, nzb_s_inner is used below, although strictly
[1]	460	!-- ---- speaking the following k-loop would have to be split up and
	461	!-- rearranged according to the staggered grid.
[132]	462	!-- However, this implies no error since staggered velocity components
	463	!-- are zero at the walls and inside buildings.
[1]	464	tn = 0
[82]	465	#if defined( __intel_openmp_bug )
[1]	466	!$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, &
	467	!$OMP tn, ust, ust2, u2, vst, vst2, v2, w2 )
	468	tn = omp_get_thread_num()
	469	#else
	470	!$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2, w2 )
	471	!$ tn = omp_get_thread_num()
	472	#endif
	473	!$OMP DO
	474	DO i = nxl, nxr
	475	DO j = nys, nyn
	476	sums_l_etot = 0.0
[132]	477	DO k = nzb_s_inner(j,i), nzt+1
[1]	478	!
	479	!-- Prognostic and diagnostic variables
	480	sums_l(k,3,tn) = sums_l(k,3,tn) + w(k,j,i) * rmask(j,i,sr)
	481	sums_l(k,8,tn) = sums_l(k,8,tn) + e(k,j,i) * rmask(j,i,sr)
	482	sums_l(k,9,tn) = sums_l(k,9,tn) + km(k,j,i) * rmask(j,i,sr)
	483	sums_l(k,10,tn) = sums_l(k,10,tn) + kh(k,j,i) * rmask(j,i,sr)
	484	sums_l(k,40,tn) = sums_l(k,40,tn) + p(k,j,i)
	485
	486	sums_l(k,33,tn) = sums_l(k,33,tn) + &
	487	( pt(k,j,i)-hom(k,1,4,sr) )*2 rmask(j,i,sr)
[624]	488
	489	IF ( humidity ) THEN
	490	sums_l(k,70,tn) = sums_l(k,70,tn) + &
	491	( q(k,j,i)-hom(k,1,41,sr) )*2 rmask(j,i,sr)
	492	ENDIF
[1007]	493
[699]	494	!
	495	!-- Higher moments
	496	!-- (Computation of the skewness of w further below)
	497	sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i)*3 rmask(j,i,sr)
[667]	498
[1]	499	sums_l_etot = sums_l_etot + &
[667]	500	0.5 * ( u(k,j,i)2 + v(k,j,i)2 + &
	501	w(k,j,i)*2 ) rmask(j,i,sr)
[1]	502	ENDDO
	503	!
	504	!-- Total and perturbation energy for the total domain (being
	505	!-- collected in the last column of sums_l). Summation of these
	506	!-- quantities is seperated from the previous loop in order to
	507	!-- allow vectorization of that loop.
[87]	508	sums_l(nzb+4,pr_palm,tn) = sums_l(nzb+4,pr_palm,tn) + sums_l_etot
[1]	509	!
	510	!-- 2D-arrays (being collected in the last column of sums_l)
[87]	511	sums_l(nzb,pr_palm,tn) = sums_l(nzb,pr_palm,tn) + &
[1]	512	us(j,i) * rmask(j,i,sr)
[87]	513	sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) + &
[1]	514	usws(j,i) * rmask(j,i,sr)
[87]	515	sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) + &
[1]	516	vsws(j,i) * rmask(j,i,sr)
[87]	517	sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) + &
[1]	518	ts(j,i) * rmask(j,i,sr)
[197]	519	IF ( humidity ) THEN
	520	sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) + &
	521	qs(j,i) * rmask(j,i,sr)
	522	ENDIF
[1]	523	ENDDO
	524	ENDDO
	525
	526	!
[667]	527	!-- Computation of statistics when ws-scheme is not used. Else these
	528	!-- quantities are evaluated in the advection routines.
[743]	529	IF ( .NOT. ws_scheme_mom .OR. sr /= 0 ) THEN
[667]	530	!$OMP DO
	531	DO i = nxl, nxr
	532	DO j = nys, nyn
	533	sums_l_eper = 0.0
	534	DO k = nzb_s_inner(j,i), nzt+1
	535	u2 = u(k,j,i)**2
	536	v2 = v(k,j,i)**2
	537	w2 = w(k,j,i)**2
	538	ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
	539	vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
	540
	541	sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr)
	542	sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr)
	543	sums_l(k,32,tn) = sums_l(k,32,tn) + w2 * rmask(j,i,sr)
	544	!
	545	!-- Perturbation energy
	546
	547	sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5 * &
	548	( ust2 + vst2 + w2 ) * rmask(j,i,sr)
	549	sums_l_eper = sums_l_eper + &
	550	0.5 * ( ust2+vst2+w2 ) * rmask(j,i,sr)
	551
	552	ENDDO
	553	sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn) &
	554	+ sums_l_eper
	555	ENDDO
	556	ENDDO
	557	ENDIF
	558	!
[1]	559	!-- Horizontally averaged profiles of the vertical fluxes
[667]	560
[1]	561	!$OMP DO
	562	DO i = nxl, nxr
	563	DO j = nys, nyn
	564	!
	565	!-- Subgridscale fluxes (without Prandtl layer from k=nzb,
	566	!-- oterwise from k=nzb+1)
[132]	567	!-- NOTE: for simplicity, nzb_diff_s_inner is used below, although
[1]	568	!-- ---- strictly speaking the following k-loop would have to be
	569	!-- split up according to the staggered grid.
[132]	570	!-- However, this implies no error since staggered velocity
	571	!-- components are zero at the walls and inside buildings.
	572
	573	DO k = nzb_diff_s_inner(j,i)-1, nzt_diff
[1]	574	!
	575	!-- Momentum flux w"u"
	576	sums_l(k,12,tn) = sums_l(k,12,tn) - 0.25 * ( &
	577	km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
	578	) * ( &
	579	( u(k+1,j,i) - u(k,j,i) ) * ddzu(k+1) &
	580	+ ( w(k,j,i) - w(k,j,i-1) ) * ddx &
	581	) * rmask(j,i,sr)
	582	!
	583	!-- Momentum flux w"v"
	584	sums_l(k,14,tn) = sums_l(k,14,tn) - 0.25 * ( &
	585	km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
	586	) * ( &
	587	( v(k+1,j,i) - v(k,j,i) ) * ddzu(k+1) &
	588	+ ( w(k,j,i) - w(k,j-1,i) ) * ddy &
	589	) * rmask(j,i,sr)
	590	!
	591	!-- Heat flux w"pt"
	592	sums_l(k,16,tn) = sums_l(k,16,tn) &
	593	- 0.5 * ( kh(k,j,i) + kh(k+1,j,i) ) &
	594	* ( pt(k+1,j,i) - pt(k,j,i) ) &
	595	* ddzu(k+1) * rmask(j,i,sr)
	596
	597
	598	!
[96]	599	!-- Salinity flux w"sa"
	600	IF ( ocean ) THEN
	601	sums_l(k,65,tn) = sums_l(k,65,tn) &
	602	- 0.5 * ( kh(k,j,i) + kh(k+1,j,i) ) &
	603	* ( sa(k+1,j,i) - sa(k,j,i) ) &
	604	* ddzu(k+1) * rmask(j,i,sr)
	605	ENDIF
	606
	607	!
[1]	608	!-- Buoyancy flux, water flux (humidity flux) w"q"
[75]	609	IF ( humidity ) THEN
[1]	610	sums_l(k,45,tn) = sums_l(k,45,tn) &
	611	- 0.5 * ( kh(k,j,i) + kh(k+1,j,i) ) &
	612	* ( vpt(k+1,j,i) - vpt(k,j,i) ) &
	613	* ddzu(k+1) * rmask(j,i,sr)
	614	sums_l(k,48,tn) = sums_l(k,48,tn) &
	615	- 0.5 * ( kh(k,j,i) + kh(k+1,j,i) ) &
	616	* ( q(k+1,j,i) - q(k,j,i) ) &
	617	* ddzu(k+1) * rmask(j,i,sr)
[1007]	618
[1]	619	IF ( cloud_physics ) THEN
	620	sums_l(k,51,tn) = sums_l(k,51,tn) &
	621	- 0.5 * ( kh(k,j,i) + kh(k+1,j,i) ) &
	622	* ( ( q(k+1,j,i) - ql(k+1,j,i) )&
	623	- ( q(k,j,i) - ql(k,j,i) ) ) &
	624	* ddzu(k+1) * rmask(j,i,sr)
	625	ENDIF
	626	ENDIF
	627
	628	!
	629	!-- Passive scalar flux
	630	IF ( passive_scalar ) THEN
	631	sums_l(k,48,tn) = sums_l(k,48,tn) &
	632	- 0.5 * ( kh(k,j,i) + kh(k+1,j,i) ) &
	633	* ( q(k+1,j,i) - q(k,j,i) ) &
	634	* ddzu(k+1) * rmask(j,i,sr)
	635	ENDIF
	636
	637	ENDDO
	638
	639	!
	640	!-- Subgridscale fluxes in the Prandtl layer
	641	IF ( use_surface_fluxes ) THEN
	642	sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
	643	usws(j,i) * rmask(j,i,sr) ! w"u"
	644	sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
	645	vsws(j,i) * rmask(j,i,sr) ! w"v"
	646	sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
	647	shf(j,i) * rmask(j,i,sr) ! w"pt"
	648	sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
	649	0.0 * rmask(j,i,sr) ! u"pt"
	650	sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
	651	0.0 * rmask(j,i,sr) ! v"pt"
[96]	652	IF ( ocean ) THEN
	653	sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
	654	saswsb(j,i) * rmask(j,i,sr) ! w"sa"
	655	ENDIF
[75]	656	IF ( humidity ) THEN
[1]	657	sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
	658	qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
[1007]	659	sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + ( &
	660	( 1.0 + 0.61 * q(nzb,j,i) ) * &
	661	shf(j,i) + 0.61 * pt(nzb,j,i) * &
	662	qsws(j,i) )
	663	IF ( cloud_droplets ) THEN
	664	sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + ( &
	665	( 1.0 + 0.61 * q(nzb,j,i) - &
	666	ql(nzb,j,i) ) * shf(j,i) + &
	667	0.61 * pt(nzb,j,i) * qsws(j,i) )
	668	ENDIF
[1]	669	IF ( cloud_physics ) THEN
	670	!
	671	!-- Formula does not work if ql(nzb) /= 0.0
	672	sums_l(nzb,51,tn) = sums_l(nzb,51,tn) + & ! w"q" (w"qv")
	673	qsws(j,i) * rmask(j,i,sr)
	674	ENDIF
	675	ENDIF
	676	IF ( passive_scalar ) THEN
	677	sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
	678	qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
	679	ENDIF
	680	ENDIF
	681
	682	!
[19]	683	!-- Subgridscale fluxes at the top surface
	684	IF ( use_top_fluxes ) THEN
[550]	685	sums_l(nzt:nzt+1,12,tn) = sums_l(nzt:nzt+1,12,tn) + &
[102]	686	uswst(j,i) * rmask(j,i,sr) ! w"u"
[550]	687	sums_l(nzt:nzt+1,14,tn) = sums_l(nzt:nzt+1,14,tn) + &
[102]	688	vswst(j,i) * rmask(j,i,sr) ! w"v"
[550]	689	sums_l(nzt:nzt+1,16,tn) = sums_l(nzt:nzt+1,16,tn) + &
[19]	690	tswst(j,i) * rmask(j,i,sr) ! w"pt"
[550]	691	sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + &
[19]	692	0.0 * rmask(j,i,sr) ! u"pt"
[550]	693	sums_l(nzt:nzt+1,61,tn) = sums_l(nzt:nzt+1,61,tn) + &
	694	0.0 * rmask(j,i,sr) ! v"pt"
	695
[96]	696	IF ( ocean ) THEN
	697	sums_l(nzt,65,tn) = sums_l(nzt,65,tn) + &
	698	saswst(j,i) * rmask(j,i,sr) ! w"sa"
	699	ENDIF
[75]	700	IF ( humidity ) THEN
[19]	701	sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
[388]	702	qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
[1007]	703	sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + ( &
	704	( 1.0 + 0.61 * q(nzt,j,i) ) * &
	705	tswst(j,i) + 0.61 * pt(nzt,j,i) * &
	706	qswst(j,i) )
	707	IF ( cloud_droplets ) THEN
	708	sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + ( &
	709	( 1.0 + 0.61 * q(nzt,j,i) - &
	710	ql(nzt,j,i) ) * tswst(j,i) + &
	711	0.61 * pt(nzt,j,i) * qswst(j,i) )
	712	ENDIF
[19]	713	IF ( cloud_physics ) THEN
	714	!
	715	!-- Formula does not work if ql(nzb) /= 0.0
	716	sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + & ! w"q" (w"qv")
	717	qswst(j,i) * rmask(j,i,sr)
	718	ENDIF
	719	ENDIF
	720	IF ( passive_scalar ) THEN
	721	sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
[388]	722	qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
[19]	723	ENDIF
	724	ENDIF
	725
	726	!
[1]	727	!-- Resolved fluxes (can be computed for all horizontal points)
[132]	728	!-- NOTE: for simplicity, nzb_s_inner is used below, although strictly
[1]	729	!-- ---- speaking the following k-loop would have to be split up and
	730	!-- rearranged according to the staggered grid.
[132]	731	DO k = nzb_s_inner(j,i), nzt
[1]	732	ust = 0.5 * ( u(k,j,i) - hom(k,1,1,sr) + &
	733	u(k+1,j,i) - hom(k+1,1,1,sr) )
	734	vst = 0.5 * ( v(k,j,i) - hom(k,1,2,sr) + &
	735	v(k+1,j,i) - hom(k+1,1,2,sr) )
	736	pts = 0.5 * ( pt(k,j,i) - hom(k,1,4,sr) + &
	737	pt(k+1,j,i) - hom(k+1,1,4,sr) )
[667]	738
[1]	739	!-- Higher moments
	740	sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)*2 &
	741	rmask(j,i,sr)
	742	sums_l(k,36,tn) = sums_l(k,36,tn) + pts*2 w(k,j,i) * &
	743	rmask(j,i,sr)
	744
	745	!
[96]	746	!-- Salinity flux and density (density does not belong to here,
[97]	747	!-- but so far there is no other suitable place to calculate)
[96]	748	IF ( ocean ) THEN
[743]	749	IF( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN
[667]	750	pts = 0.5 * ( sa(k,j,i) - hom(k,1,23,sr) + &
[96]	751	sa(k+1,j,i) - hom(k+1,1,23,sr) )
[667]	752	sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) * &
[96]	753	rmask(j,i,sr)
[667]	754	ENDIF
[96]	755	sums_l(k,64,tn) = sums_l(k,64,tn) + rho(k,j,i) * &
	756	rmask(j,i,sr)
[388]	757	sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) * &
	758	rmask(j,i,sr)
[96]	759	ENDIF
	760
	761	!
[1053]	762	!-- Buoyancy flux, water flux, humidity flux, liquid water
	763	!-- content, rain drop concentration and rain water content
[75]	764	IF ( humidity ) THEN
[1007]	765	IF ( cloud_physics .OR. cloud_droplets ) THEN
[1053]	766	pts = 0.5 * ( vpt(k,j,i) - hom(k,1,44,sr) + &
[1007]	767	vpt(k+1,j,i) - hom(k+1,1,44,sr) )
	768	sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
[1]	769	rmask(j,i,sr)
[1053]	770	IF ( .NOT. cloud_droplets ) THEN
[1115]	771	pts = 0.5 * &
	772	( ( q(k,j,i) - ql(k,j,i) ) - &
	773	hom(k,1,42,sr) + &
	774	( q(k+1,j,i) - ql(k+1,j,i) ) - &
[1053]	775	hom(k+1,1,42,sr) )
[1115]	776	sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) * &
[1053]	777	rmask(j,i,sr)
	778	IF ( icloud_scheme == 0 ) THEN
[1115]	779	sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * &
[1053]	780	rmask(j,i,sr)
[1115]	781	sums_l(k,75,tn) = sums_l(k,75,tn) + qc(k,j,i) * &
[1053]	782	rmask(j,i,sr)
[1115]	783	IF ( precipitation ) THEN
	784	sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) * &
	785	rmask(j,i,sr)
	786	sums_l(k,74,tn) = sums_l(k,74,tn) + qr(k,j,i) * &
	787	rmask(j,i,sr)
	788	sums_l(k,76,tn) = sums_l(k,76,tn) + prr(k,j,i) *&
	789	rmask(j,i,sr)
	790	ENDIF
[1053]	791	ELSE
[1115]	792	sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * &
[1053]	793	rmask(j,i,sr)
	794	ENDIF
	795	ELSE
[1115]	796	sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * &
[1053]	797	rmask(j,i,sr)
	798	ENDIF
[1007]	799	ELSE
	800	IF( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN
	801	pts = 0.5 * ( vpt(k,j,i) - hom(k,1,44,sr) + &
	802	vpt(k+1,j,i) - hom(k+1,1,44,sr) )
	803	sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
	804	rmask(j,i,sr)
	805	ELSE IF ( ws_scheme_sca .AND. sr == 0 ) THEN
[1053]	806	sums_l(k,46,tn) = ( 1.0 + 0.61 * hom(k,1,41,sr) ) * &
[1007]	807	sums_l(k,17,tn) + &
	808	0.61 * hom(k,1,4,sr) * sums_l(k,49,tn)
	809	END IF
	810	END IF
[1]	811	ENDIF
	812	!
	813	!-- Passive scalar flux
[743]	814	IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca &
	815	.OR. sr /= 0 ) ) THEN
[1]	816	pts = 0.5 * ( q(k,j,i) - hom(k,1,41,sr) + &
	817	q(k+1,j,i) - hom(k+1,1,41,sr) )
	818	sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
	819	rmask(j,i,sr)
	820	ENDIF
	821
	822	!
	823	!-- Energy flux we
[667]	824	!-- has to be adjusted
	825	sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5 * &
	826	( ust2 + vst2 + w(k,j,i)**2 )&
	827	* rmask(j,i,sr)
[1]	828	ENDDO
	829	ENDDO
	830	ENDDO
[709]	831	!
	832	!-- For speed optimization fluxes which have been computed in part directly
	833	!-- inside the WS advection routines are treated seperatly
	834	!-- Momentum fluxes first:
[743]	835	IF ( .NOT. ws_scheme_mom .OR. sr /= 0 ) THEN
[667]	836	!$OMP DO
	837	DO i = nxl, nxr
	838	DO j = nys, nyn
	839	DO k = nzb_diff_s_inner(j,i)-1, nzt_diff
	840	ust = 0.5 * ( u(k,j,i) - hom(k,1,1,sr) + &
	841	u(k+1,j,i) - hom(k+1,1,1,sr) )
	842	vst = 0.5 * ( v(k,j,i) - hom(k,1,2,sr) + &
	843	v(k+1,j,i) - hom(k+1,1,2,sr) )
[1007]	844	!
[667]	845	!-- Momentum flux wu
	846	sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5 * &
	847	( w(k,j,i-1) + w(k,j,i) ) &
	848	* ust * rmask(j,i,sr)
	849	!
	850	!-- Momentum flux wv
	851	sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5 * &
	852	( w(k,j-1,i) + w(k,j,i) ) &
	853	* vst * rmask(j,i,sr)
	854	ENDDO
	855	ENDDO
	856	ENDDO
[1]	857
[667]	858	ENDIF
[743]	859	IF ( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN
[667]	860	!$OMP DO
	861	DO i = nxl, nxr
	862	DO j = nys, nyn
[709]	863	DO k = nzb_diff_s_inner(j,i)-1, nzt_diff
	864	!
	865	!-- Vertical heat flux
[667]	866	sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5 * &
	867	( pt(k,j,i) - hom(k,1,4,sr) + &
	868	pt(k+1,j,i) - hom(k+1,1,4,sr) ) &
	869	* w(k,j,i) * rmask(j,i,sr)
	870	IF ( humidity ) THEN
	871	pts = 0.5 * ( q(k,j,i) - hom(k,1,41,sr) + &
	872	q(k+1,j,i) - hom(k+1,1,41,sr) )
	873	sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
	874	rmask(j,i,sr)
	875	ENDIF
	876	ENDDO
	877	ENDDO
	878	ENDDO
	879
	880	ENDIF
	881
	882
[1]	883	!
[97]	884	!-- Density at top follows Neumann condition
[388]	885	IF ( ocean ) THEN
	886	sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
	887	sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
	888	ENDIF
[97]	889
	890	!
[1]	891	!-- Divergence of vertical flux of resolved scale energy and pressure
[106]	892	!-- fluctuations as well as flux of pressure fluctuation itself (68).
	893	!-- First calculate the products, then the divergence.
[1]	894	!-- Calculation is time consuming. Do it only, if profiles shall be plotted.
[106]	895	IF ( hom(nzb+1,2,55,0) /= 0.0 .OR. hom(nzb+1,2,68,0) /= 0.0 ) THEN
[1]	896
	897	sums_ll = 0.0 ! local array
	898
	899	!$OMP DO
	900	DO i = nxl, nxr
	901	DO j = nys, nyn
[132]	902	DO k = nzb_s_inner(j,i)+1, nzt
[1]	903
	904	sums_ll(k,1) = sums_ll(k,1) + 0.5 * w(k,j,i) * ( &
	905	( 0.25 * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1) &
[678]	906	- 0.5 * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) ) &
[1]	907	) )**2 &
	908	+ ( 0.25 * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i) &
[678]	909	- 0.5 * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) ) &
[1]	910	) )**2 &
	911	+ w(k,j,i)**2 )
	912
	913	sums_ll(k,2) = sums_ll(k,2) + 0.5 * w(k,j,i) &
	914	* ( p(k,j,i) + p(k+1,j,i) )
	915
	916	ENDDO
	917	ENDDO
	918	ENDDO
	919	sums_ll(0,1) = 0.0 ! because w is zero at the bottom
	920	sums_ll(nzt+1,1) = 0.0
	921	sums_ll(0,2) = 0.0
	922	sums_ll(nzt+1,2) = 0.0
	923
[678]	924	DO k = nzb+1, nzt
[1]	925	sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
	926	sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
[106]	927	sums_l(k,68,tn) = sums_ll(k,2)
[1]	928	ENDDO
	929	sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
	930	sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
[106]	931	sums_l(nzb,68,tn) = 0.0 ! because w* = 0 at nzb
[1]	932
	933	ENDIF
	934
	935	!
[106]	936	!-- Divergence of vertical flux of SGS TKE and the flux itself (69)
	937	IF ( hom(nzb+1,2,57,0) /= 0.0 .OR. hom(nzb+1,2,69,0) /= 0.0 ) THEN
[1]	938
	939	!$OMP DO
	940	DO i = nxl, nxr
	941	DO j = nys, nyn
[132]	942	DO k = nzb_s_inner(j,i)+1, nzt
[1]	943
[106]	944	sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5 * ( &
[1]	945	(km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
	946	- (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k) &
[106]	947	) * ddzw(k)
[1]	948
[106]	949	sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5 * ( &
	950	(km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
	951	)
	952
[1]	953	ENDDO
	954	ENDDO
	955	ENDDO
	956	sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
[106]	957	sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
[1]	958
	959	ENDIF
	960
	961	!
	962	!-- Horizontal heat fluxes (subgrid, resolved, total).
	963	!-- Do it only, if profiles shall be plotted.
	964	IF ( hom(nzb+1,2,58,0) /= 0.0 ) THEN
	965
	966	!$OMP DO
	967	DO i = nxl, nxr
	968	DO j = nys, nyn
[132]	969	DO k = nzb_s_inner(j,i)+1, nzt
[1]	970	!
	971	!-- Subgrid horizontal heat fluxes u"pt", v"pt"
	972	sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5 * &
	973	( kh(k,j,i) + kh(k,j,i-1) ) &
	974	* ( pt(k,j,i-1) - pt(k,j,i) ) &
	975	* ddx * rmask(j,i,sr)
	976	sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5 * &
	977	( kh(k,j,i) + kh(k,j-1,i) ) &
	978	* ( pt(k,j-1,i) - pt(k,j,i) ) &
	979	* ddy * rmask(j,i,sr)
	980	!
	981	!-- Resolved horizontal heat fluxes upt, vpt
	982	sums_l(k,59,tn) = sums_l(k,59,tn) + &
	983	( u(k,j,i) - hom(k,1,1,sr) ) &
	984	* 0.5 * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
	985	pt(k,j,i) - hom(k,1,4,sr) )
	986	pts = 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
	987	pt(k,j,i) - hom(k,1,4,sr) )
	988	sums_l(k,62,tn) = sums_l(k,62,tn) + &
	989	( v(k,j,i) - hom(k,1,2,sr) ) &
	990	* 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
	991	pt(k,j,i) - hom(k,1,4,sr) )
	992	ENDDO
	993	ENDDO
	994	ENDDO
	995	!
	996	!-- Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
[97]	997	sums_l(nzb,58,tn) = 0.0
	998	sums_l(nzb,59,tn) = 0.0
	999	sums_l(nzb,60,tn) = 0.0
	1000	sums_l(nzb,61,tn) = 0.0
	1001	sums_l(nzb,62,tn) = 0.0
	1002	sums_l(nzb,63,tn) = 0.0
[1]	1003
	1004	ENDIF
[87]	1005
	1006	!
	1007	!-- Calculate the user-defined profiles
	1008	CALL user_statistics( 'profiles', sr, tn )
[1]	1009	!$OMP END PARALLEL
	1010
	1011	!
	1012	!-- Summation of thread sums
	1013	IF ( threads_per_task > 1 ) THEN
	1014	DO i = 1, threads_per_task-1
	1015	sums_l(:,3,0) = sums_l(:,3,0) + sums_l(:,3,i)
	1016	sums_l(:,4:40,0) = sums_l(:,4:40,0) + sums_l(:,4:40,i)
[87]	1017	sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
	1018	sums_l(:,45:pr_palm,i)
	1019	IF ( max_pr_user > 0 ) THEN
	1020	sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
	1021	sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
	1022	sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
	1023	ENDIF
[1]	1024	ENDDO
	1025	ENDIF
	1026
	1027	#if defined( __parallel )
[667]	1028
[1]	1029	!
	1030	!-- Compute total sum from local sums
[622]	1031	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[1]	1032	CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, &
	1033	MPI_SUM, comm2d, ierr )
	1034	#else
	1035	sums = sums_l(:,:,0)
	1036	#endif
	1037
	1038	!
	1039	!-- Final values are obtained by division by the total number of grid points
	1040	!-- used for summation. After that store profiles.
	1041	!-- Profiles:
	1042	DO k = nzb, nzt+1
[132]	1043	sums(k,3) = sums(k,3) / ngp_2dh(sr)
[142]	1044	sums(k,8:11) = sums(k,8:11) / ngp_2dh_s_inner(k,sr)
[132]	1045	sums(k,12:22) = sums(k,12:22) / ngp_2dh(sr)
	1046	sums(k,23:29) = sums(k,23:29) / ngp_2dh_s_inner(k,sr)
	1047	sums(k,30:32) = sums(k,30:32) / ngp_2dh(sr)
[142]	1048	sums(k,33:34) = sums(k,33:34) / ngp_2dh_s_inner(k,sr)
	1049	sums(k,35:39) = sums(k,35:39) / ngp_2dh(sr)
[132]	1050	sums(k,40) = sums(k,40) / ngp_2dh_s_inner(k,sr)
	1051	sums(k,45:53) = sums(k,45:53) / ngp_2dh(sr)
	1052	sums(k,54) = sums(k,54) / ngp_2dh_s_inner(k,sr)
	1053	sums(k,55:63) = sums(k,55:63) / ngp_2dh(sr)
	1054	sums(k,64) = sums(k,64) / ngp_2dh_s_inner(k,sr)
	1055	sums(k,65:69) = sums(k,65:69) / ngp_2dh(sr)
	1056	sums(k,70:pr_palm-2) = sums(k,70:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
[1]	1057	ENDDO
[667]	1058
[1]	1059	!-- Upstream-parts
[87]	1060	sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
[1]	1061	!-- u* and so on
[87]	1062	!-- As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
[1]	1063	!-- size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
	1064	!-- above the topography, they are being divided by ngp_2dh(sr)
[87]	1065	sums(nzb:nzb+3,pr_palm) = sums(nzb:nzb+3,pr_palm) / &
[1]	1066	ngp_2dh(sr)
[197]	1067	sums(nzb+12,pr_palm) = sums(nzb+12,pr_palm) / & ! qs
	1068	ngp_2dh(sr)
[1]	1069	!-- eges, e*
[87]	1070	sums(nzb+4:nzb+5,pr_palm) = sums(nzb+4:nzb+5,pr_palm) / &
[132]	1071	ngp_3d(sr)
[1]	1072	!-- Old and new divergence
[87]	1073	sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
[1]	1074	ngp_3d_inner(sr)
	1075
[87]	1076	!-- User-defined profiles
	1077	IF ( max_pr_user > 0 ) THEN
	1078	DO k = nzb, nzt+1
	1079	sums(k,pr_palm+1:pr_palm+max_pr_user) = &
	1080	sums(k,pr_palm+1:pr_palm+max_pr_user) / &
[132]	1081	ngp_2dh_s_inner(k,sr)
[87]	1082	ENDDO
	1083	ENDIF
[1007]	1084
[1]	1085	!
	1086	!-- Collect horizontal average in hom.
	1087	!-- Compute deduced averages (e.g. total heat flux)
	1088	hom(:,1,3,sr) = sums(:,3) ! w
	1089	hom(:,1,8,sr) = sums(:,8) ! e profiles 5-7 are initial profiles
	1090	hom(:,1,9,sr) = sums(:,9) ! km
	1091	hom(:,1,10,sr) = sums(:,10) ! kh
	1092	hom(:,1,11,sr) = sums(:,11) ! l
	1093	hom(:,1,12,sr) = sums(:,12) ! w"u"
	1094	hom(:,1,13,sr) = sums(:,13) ! wu
	1095	hom(:,1,14,sr) = sums(:,14) ! w"v"
	1096	hom(:,1,15,sr) = sums(:,15) ! wv
	1097	hom(:,1,16,sr) = sums(:,16) ! w"pt"
	1098	hom(:,1,17,sr) = sums(:,17) ! wpt
	1099	hom(:,1,18,sr) = sums(:,16) + sums(:,17) ! wpt
	1100	hom(:,1,19,sr) = sums(:,12) + sums(:,13) ! wu
	1101	hom(:,1,20,sr) = sums(:,14) + sums(:,15) ! wv
	1102	hom(:,1,21,sr) = sums(:,21) ! wptBC
	1103	hom(:,1,22,sr) = sums(:,16) + sums(:,21) ! wptBC
[96]	1104	! profile 24 is initial profile (sa)
	1105	! profiles 25-29 left empty for initial
[1]	1106	! profiles
	1107	hom(:,1,30,sr) = sums(:,30) ! u*2
	1108	hom(:,1,31,sr) = sums(:,31) ! v*2
	1109	hom(:,1,32,sr) = sums(:,32) ! w*2
	1110	hom(:,1,33,sr) = sums(:,33) ! pt*2
	1111	hom(:,1,34,sr) = sums(:,34) ! e*
	1112	hom(:,1,35,sr) = sums(:,35) ! w2pt
	1113	hom(:,1,36,sr) = sums(:,36) ! wpt2
	1114	hom(:,1,37,sr) = sums(:,37) ! we
	1115	hom(:,1,38,sr) = sums(:,38) ! w*3
[699]	1116	hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20 )**1.5 ! Sw
[1]	1117	hom(:,1,40,sr) = sums(:,40) ! p
[531]	1118	hom(:,1,45,sr) = sums(:,45) ! w"vpt"
[1]	1119	hom(:,1,46,sr) = sums(:,46) ! wvpt
	1120	hom(:,1,47,sr) = sums(:,45) + sums(:,46) ! wvpt
	1121	hom(:,1,48,sr) = sums(:,48) ! w"q" (w"qv")
	1122	hom(:,1,49,sr) = sums(:,49) ! wq (wqv)
	1123	hom(:,1,50,sr) = sums(:,48) + sums(:,49) ! wq (wqv)
	1124	hom(:,1,51,sr) = sums(:,51) ! w"qv"
	1125	hom(:,1,52,sr) = sums(:,52) ! wqv
	1126	hom(:,1,53,sr) = sums(:,52) + sums(:,51) ! wq (wqv)
	1127	hom(:,1,54,sr) = sums(:,54) ! ql
	1128	hom(:,1,55,sr) = sums(:,55) ! wuu*/dz
	1129	hom(:,1,56,sr) = sums(:,56) ! wp/dz
[106]	1130	hom(:,1,57,sr) = sums(:,57) ! ( w"e + w"p"/rho )/dz
[1]	1131	hom(:,1,58,sr) = sums(:,58) ! u"pt"
	1132	hom(:,1,59,sr) = sums(:,59) ! upt
	1133	hom(:,1,60,sr) = sums(:,58) + sums(:,59) ! upt_t
	1134	hom(:,1,61,sr) = sums(:,61) ! v"pt"
	1135	hom(:,1,62,sr) = sums(:,62) ! vpt
	1136	hom(:,1,63,sr) = sums(:,61) + sums(:,62) ! vpt_t
[96]	1137	hom(:,1,64,sr) = sums(:,64) ! rho
	1138	hom(:,1,65,sr) = sums(:,65) ! w"sa"
	1139	hom(:,1,66,sr) = sums(:,66) ! wsa
	1140	hom(:,1,67,sr) = sums(:,65) + sums(:,66) ! wsa
[106]	1141	hom(:,1,68,sr) = sums(:,68) ! wp
	1142	hom(:,1,69,sr) = sums(:,69) ! w"e + w"p"/rho
[197]	1143	hom(:,1,70,sr) = sums(:,70) ! q*2
[388]	1144	hom(:,1,71,sr) = sums(:,71) ! prho
[531]	1145	hom(:,1,72,sr) = hyp * 1E-4 ! hyp in dbar
[1053]	1146	hom(:,1,73,sr) = sums(:,73) ! nr
	1147	hom(:,1,74,sr) = sums(:,74) ! qr
	1148	hom(:,1,75,sr) = sums(:,75) ! qc
	1149	hom(:,1,76,sr) = sums(:,76) ! prr (precipitation rate)
[1179]	1150	! 77 is initial density profile
[1]	1151
[87]	1152	hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
[1]	1153	! upstream-parts u_x, u_y, u_z, v_x,
	1154	! v_y, usw. (in last but one profile)
[667]	1155	hom(:,1,pr_palm,sr) = sums(:,pr_palm)
[1]	1156	! u, w'u', w'v', t (in last profile)
	1157
[87]	1158	IF ( max_pr_user > 0 ) THEN ! user-defined profiles
	1159	hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
	1160	sums(:,pr_palm+1:pr_palm+max_pr_user)
	1161	ENDIF
	1162
[1]	1163	!
	1164	!-- Determine the boundary layer height using two different schemes.
[94]	1165	!-- First scheme: Starting from the Earth's (Ocean's) surface, look for the
	1166	!-- first relative minimum (maximum) of the total heat flux.
	1167	!-- The corresponding height is assumed as the boundary layer height, if it
	1168	!-- is less than 1.5 times the height where the heat flux becomes negative
	1169	!-- (positive) for the first time.
[1]	1170	z_i(1) = 0.0
	1171	first = .TRUE.
[667]	1172
[97]	1173	IF ( ocean ) THEN
	1174	DO k = nzt, nzb+1, -1
[667]	1175	IF ( first .AND. hom(k,1,18,sr) < 0.0 &
	1176	.AND. abs(hom(k,1,18,sr)) > 1.0E-8) THEN
[97]	1177	first = .FALSE.
	1178	height = zw(k)
	1179	ENDIF
	1180	IF ( hom(k,1,18,sr) < 0.0 .AND. &
[667]	1181	abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
[97]	1182	hom(k-1,1,18,sr) > hom(k,1,18,sr) ) THEN
	1183	IF ( zw(k) < 1.5 * height ) THEN
	1184	z_i(1) = zw(k)
	1185	ELSE
	1186	z_i(1) = height
	1187	ENDIF
	1188	EXIT
	1189	ENDIF
	1190	ENDDO
	1191	ELSE
[94]	1192	DO k = nzb, nzt-1
[667]	1193	IF ( first .AND. hom(k,1,18,sr) < 0.0 &
	1194	.AND. abs(hom(k,1,18,sr)) > 1.0E-8 ) THEN
[94]	1195	first = .FALSE.
	1196	height = zw(k)
[1]	1197	ENDIF
[667]	1198	IF ( hom(k,1,18,sr) < 0.0 .AND. &
	1199	abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
[94]	1200	hom(k+1,1,18,sr) > hom(k,1,18,sr) ) THEN
	1201	IF ( zw(k) < 1.5 * height ) THEN
	1202	z_i(1) = zw(k)
	1203	ELSE
	1204	z_i(1) = height
	1205	ENDIF
	1206	EXIT
	1207	ENDIF
	1208	ENDDO
[97]	1209	ENDIF
[1]	1210
	1211	!
[291]	1212	!-- Second scheme: Gradient scheme from Sullivan et al. (1998), modified
	1213	!-- by Uhlenbrock(2006). The boundary layer height is the height with the
	1214	!-- maximal local temperature gradient: starting from the second (the last
	1215	!-- but one) vertical gridpoint, the local gradient must be at least
	1216	!-- 0.2K/100m and greater than the next four gradients.
	1217	!-- WARNING: The threshold value of 0.2K/100m must be adjusted for the
	1218	!-- ocean case!
[1]	1219	z_i(2) = 0.0
[291]	1220	DO k = nzb+1, nzt+1
	1221	dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
	1222	ENDDO
	1223	dptdz_threshold = 0.2 / 100.0
	1224
[97]	1225	IF ( ocean ) THEN
[291]	1226	DO k = nzt+1, nzb+5, -1
	1227	IF ( dptdz(k) > dptdz_threshold .AND. &
	1228	dptdz(k) > dptdz(k-1) .AND. dptdz(k) > dptdz(k-2) .AND. &
	1229	dptdz(k) > dptdz(k-3) .AND. dptdz(k) > dptdz(k-4) ) THEN
	1230	z_i(2) = zw(k-1)
[97]	1231	EXIT
	1232	ENDIF
	1233	ENDDO
	1234	ELSE
[291]	1235	DO k = nzb+1, nzt-3
	1236	IF ( dptdz(k) > dptdz_threshold .AND. &
	1237	dptdz(k) > dptdz(k+1) .AND. dptdz(k) > dptdz(k+2) .AND. &
	1238	dptdz(k) > dptdz(k+3) .AND. dptdz(k) > dptdz(k+4) ) THEN
	1239	z_i(2) = zw(k-1)
[97]	1240	EXIT
	1241	ENDIF
	1242	ENDDO
	1243	ENDIF
[1]	1244
[87]	1245	hom(nzb+6,1,pr_palm,sr) = z_i(1)
	1246	hom(nzb+7,1,pr_palm,sr) = z_i(2)
[1]	1247
	1248	!
	1249	!-- Computation of both the characteristic vertical velocity and
	1250	!-- the characteristic convective boundary layer temperature.
	1251	!-- The horizontal average at nzb+1 is input for the average temperature.
[667]	1252	IF ( hom(nzb,1,18,sr) > 0.0 .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8 &
	1253	.AND. z_i(1) /= 0.0 ) THEN
[87]	1254	hom(nzb+8,1,pr_palm,sr) = ( g / hom(nzb+1,1,4,sr) * &
[94]	1255	hom(nzb,1,18,sr) * &
	1256	ABS( z_i(1) ) )**0.333333333
[1]	1257	!-- so far this only works if Prandtl layer is used
[87]	1258	hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
[1]	1259	ELSE
[87]	1260	hom(nzb+8,1,pr_palm,sr) = 0.0
	1261	hom(nzb+11,1,pr_palm,sr) = 0.0
[1]	1262	ENDIF
	1263
[48]	1264	!
	1265	!-- Collect the time series quantities
[87]	1266	ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr) ! E
	1267	ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr) ! E*
[48]	1268	ts_value(3,sr) = dt_3d
[87]	1269	ts_value(4,sr) = hom(nzb,1,pr_palm,sr) ! u*
	1270	ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr) ! th*
[48]	1271	ts_value(6,sr) = u_max
	1272	ts_value(7,sr) = v_max
	1273	ts_value(8,sr) = w_max
[87]	1274	ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr) ! new divergence
	1275	ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr) ! old Divergence
	1276	ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr) ! z_i(1)
	1277	ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr) ! z_i(2)
	1278	ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr) ! w*
[48]	1279	ts_value(14,sr) = hom(nzb,1,16,sr) ! w'pt' at k=0
	1280	ts_value(15,sr) = hom(nzb+1,1,16,sr) ! w'pt' at k=1
	1281	ts_value(16,sr) = hom(nzb+1,1,18,sr) ! wpt at k=1
	1282	ts_value(17,sr) = hom(nzb,1,4,sr) ! pt(0)
	1283	ts_value(18,sr) = hom(nzb+1,1,4,sr) ! pt(zp)
[197]	1284	ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr) ! u'w' at k=0
	1285	ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr) ! v'w' at k=0
[343]	1286	ts_value(21,sr) = hom(nzb,1,48,sr) ! w"q" at k=0
[197]	1287
[48]	1288	IF ( ts_value(5,sr) /= 0.0 ) THEN
	1289	ts_value(22,sr) = ts_value(4,sr)**2 / &
	1290	( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
	1291	ELSE
	1292	ts_value(22,sr) = 10000.0
	1293	ENDIF
[1]	1294
[343]	1295	ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr) ! q*
[1]	1296	!
[48]	1297	!-- Calculate additional statistics provided by the user interface
[87]	1298	CALL user_statistics( 'time_series', sr, 0 )
[1]	1299
[48]	1300	ENDDO ! loop of the subregions
	1301
[1]	1302	!
	1303	!-- If required, sum up horizontal averages for subsequent time averaging
	1304	IF ( do_sum ) THEN
	1305	IF ( average_count_pr == 0 ) hom_sum = 0.0
	1306	hom_sum = hom_sum + hom(:,1,:,:)
	1307	average_count_pr = average_count_pr + 1
	1308	do_sum = .FALSE.
	1309	ENDIF
	1310
	1311	!
	1312	!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
	1313	!-- This flag is reset after each time step in time_integration.
	1314	flow_statistics_called = .TRUE.
	1315
	1316	CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
	1317
	1318
	1319	END SUBROUTINE flow_statistics

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

Context Navigation

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

Download in other formats: