Home

Context Navigation

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

Last change on this file since 681 was 679, checked in by raasch, 14 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 51.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |