Home

Context Navigation

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

Last change on this file since 76 was 75, checked in by raasch, 17 years ago
preliminary update for changes concerning non-cyclic boundary conditions
Property svn:keywords set to `Id`
File size: 36.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |