Home

Context Navigation

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

Last change on this file since 844 was 802, checked in by suehring, 13 years ago
Last commit documented.
Property svn:keywords set to `Id`
File size: 52.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |