Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |