Home

Context Navigation

flow_statistics.f90 @ 83

Last change on this file since 83 was 83, checked in by raasch, 18 years ago

New:
---

Changed:

PALM can be generally installed on any kind of Linux-, IBM-AIX-, or NEC-SX-system by adding appropriate settings to the configuration file.

Scripts are also running under the public domain ksh.

All system relevant compile and link options as well as the host identifier (local_host) are specified in the configuration file.

Filetransfer by ftp removed (options -f removed from mrun and mbuild).

Call of (system-)FLUSH routine moved to new routine local_flush.

return_addres and return_username are read from ENVPAR-NAMELIST-file instead of using local_getenv.

Preprocessor strings for different linux clusters changed to "lc", some preprocessor directives renamed (new: intel_openmp_bug), preprocessor directives for old systems removed

advec_particles, check_open, cpu_log, cpu_statistics, data_output_dvrp, flow_statistics, header, init_dvrp, init_particles, init_1d_model, init_dvrp, init_pegrid, local_getenv, local_system, local_tremain, local_tremain_ini, modules, palm, parin, run_control

new:
local_flush

mbuild, mrun

Errors:

Property svn:keywords set to Id

File size: 36.9 KB

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

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

Context Navigation

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

Download in other formats: