Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |