Home

Context Navigation

flow_statistics.f90 @ 132

Last change on this file since 132 was 132, checked in by letzel, 16 years ago

Vertical profiles now based on nzb_s_inner; they are divided by
ngp_2dh_s_inner (scalars, procucts of scalars and velocity components) and
ngp_2dh (staggered velocity components and their products), respectively.

Allow new case bc_uv_t = 'dirichlet_0' for channel flow.

Property svn:keywords set to Id

File size: 44.2 KB

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

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

Context Navigation

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

Download in other formats: