Home

Context Navigation

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

Last change on this file since 1195 was 1182, checked in by raasch, 11 years ago
last commit documented, rc-file for example run updated
Property svn:keywords set to `Id`
File size: 56.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |