Home

Context Navigation

flow_statistics.f90 @ 80

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

New:
---

particle reflection from vertical walls implemented, particle SGS model adjusted to walls

Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.

new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files

new d3par-parameter dt_max to define the maximum value for the allowed timestep

new inipar-parameter loop_optimization to control the loop optimization method

new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).

new user interface user_advec_particles

new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays

topography height informations are stored on arrays zu_s_inner and zw_w_inner and output to the 2d/3d NetCDF files

samples added to the user interface which show how to add user-define time series quantities.

calculation/output of precipitation amount, precipitation rate and z0 (by setting "pra*", "prr*", "z0*" with data_output). The time interval on which the precipitation amount is defined is set by new d3par-parameter precipitation_amount_interval

unit 9 opened for debug output (file DEBUG_<pe#>)

Makefile, advec_particles, average_3d_data, buoyancy, calc_precipitation, check_open, check_parameters, data_output_2d, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, impact_of_latent_heat, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, read_3d_binary, sum_up_3d_data, user_interface, write_var_list, write_3d_binary

New: wall_fluxes

Changed:

General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, volume flow control is added for the outflow boundary (currently only for the north boundary!!), additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d

The FFT-method for solving the Poisson-equation is now working with Neumann boundary conditions both at the bottom and the top. This requires adjustments of the tridiagonal coefficients and subtracting the horizontally averaged mean from the vertical velocity field.

+age_m in particle_type

Particles-package is now part of the default code ("-p particles" is not needed any more).

Move call of user_actions( 'after_integration' ) below increment of times
and counters. user_actions is now called for each statistic region and has as an argument the number of the respective region (sr)

d3par-parameter data_output_ts removed. Timeseries output for "profil" removed. Timeseries are now switched on by dt_dots. Timeseries data is collected in flow_statistics.

Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.

q is not allowed to become negative (prognostic_equations).

poisfft_init is only called if fft-solver is switched on (init_pegrid).

d3par-parameter moisture renamed to humidity.

Subversion global revision number is read from mrun and added to the run description header and to the run control (_rc) file.

vtk directives removed from main program.

The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.

advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, asselin_filter, check_parameters, coriolis, data_output_dvrp, data_output_ptseries, data_output_ts, data_output_2d, data_output_3d, diffusion_u, diffusion_v, exchange_horiz, exchange_horiz_2d, flow_statistics, header, init_grid, init_particles, init_pegrid, init_rankine, init_pt_anomaly, init_1d_model, init_3d_model, modules, palm, package_parin, parin, poisfft, poismg, prandtl_fluxes, pres, production_e, prognostic_equations, read_var_list, read_3d_binary, sor, swap_timelevel, time_integration, write_var_list, write_3d_binary

Errors:

Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model

Bugfix in sample for reading user defined data from restart file (user_init)

Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model)

Check for possible negative humidities in the initial humidity profile.

in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in
# case of .mod files

Makefile
check_parameters, init_1d_model, init_3d_model, user_interface

Property svn:keywords set to Id

File size: 36.8 KB

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

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

Context Navigation

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

Download in other formats: