Home

Context Navigation

flow_statistics.f90 @ 667

Last change on this file since 667 was 667, checked in by suehring, 13 years ago

summary:

Gryschka:

Coupling with different resolution and different numbers of PEs in ocean and atmosphere is available
Exchange of u and v from ocean surface to atmosphere surface
Mirror boundary condition for u and v at the bottom are replaced by dirichlet boundary conditions
Inflow turbulence is now defined by flucuations around spanwise mean
Bugfixes for cyclic_fill and constant_volume_flow

Suehring:

New advection added ( Wicker and Skamarock 5th order ), therefore:
- New module advec_ws.f90
- Modified exchange of ghost boundaries.
- Modified evaluation of turbulent fluxes
- New index bounds nxlg, nxrg, nysg, nyng

advec_ws.f90

Advection scheme for scalars and momentum using the flux formulation of
Wicker and Skamarock 5th order.
Additionally the module contains of a routine using for initialisation and
steering of the statical evaluation. The computation of turbulent fluxes takes
place inside the advection routines.
In case of vector architectures Dirichlet and Radiation boundary conditions are
outstanding and not available. Furthermore simulations within topography are
not possible so far. A further routine local_diss_ij is available and is used
if a control of dissipative fluxes is desired.

check_parameters.f90

Exchange of parameters between ocean and atmosphere via PE0
Check for illegal combination of ws-scheme and timestep scheme.
Check for topography and ws-scheme.
Check for not cyclic boundary conditions in combination with ws-scheme and
loop_optimization = 'vector'.
Check for call_psolver_at_all_substeps and ws-scheme for momentum_advec.

Different processor/grid topology in atmosphere and ocean is now allowed!
Bugfixes in checking for conserve_volume_flow_mode.

exchange_horiz.f90

Dynamic exchange of ghost points with nbgp_local to ensure that no useless
ghost points exchanged in case of multigrid. type_yz(0) and type_xz(0) used for
normal grid, the remaining types used for the several grid levels.
Exchange is done via MPI-Vectors with a dynamic value of ghost points which
depend on the advection scheme. Exchange of left and right PEs is 10% faster
with MPI-Vectors than without.

flow_statistics.f90

When advection is computed with ws-scheme, turbulent fluxes are already
computed in the respective advection routines and buffered in arrays
sums_xxxx_ws_l(). This is due to a consistent treatment of statistics
with the numerics and to avoid unphysical kinks near the surface. So some if-
requests has to be done to dicern between fluxes from ws-scheme other advection
schemes. Furthermore the computation of z_i is only done if the heat flux
exceeds a minimum value. This affects only simulations of a neutral boundary
layer and is due to reasons of computations in the advection scheme.

inflow_turbulence.f90

Using nbgp recycling planes for a better resolution of the turbulent flow near
the inflow.

init_grid.f90

Definition of new array bounds nxlg, nxrg, nysg, nyng on each PE.
Furthermore the allocation of arrays and steering of loops is done with these
parameters. Call of exchange_horiz are modified.
In case of dirichlet bounday condition at the bottom zu(0)=0.0
dzu_mg has to be set explicitly for a equally spaced grid near bottom.
ddzu_pres added to use a equally spaced grid near bottom.

init_pegrid.f90

Moved determination of target_id's from init_coupling
Determination of parameters needed for coupling (coupling_topology, ngp_a, ngp_o)
with different grid/processor-topology in ocean and atmosphere

Adaption of ngp_xy, ngp_y to a dynamic number of ghost points.
The maximum_grid_level changed from 1 to 0. 0 is the normal grid, 1 to
maximum_grid_level the grids for multigrid, in which 0 and 1 are normal grids.
This distinction is due to reasons of data exchange and performance for the
normal grid and grids in poismg.
The definition of MPI-Vectors adapted to a dynamic numer of ghost points.
New MPI-Vectors for data exchange between left and right boundaries added.
This is due to reasons of performance (10% faster).

ATTENTION: nnz_x undefined problem still has to be solved!!!!!!!!
TEST OUTPUT (TO BE REMOVED) logging mpi2 ierr values

parin.f90

Steering parameter dissipation_control added in inipar.

Makefile

Module advec_ws added.

Modules

Removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc

For coupling with different resolution in ocean and atmophere:
+nx_a, +nx_o, ny_a, +ny_o, ngp_a, ngp_o, +total_2d_o, +total_2d_a,
+coupling_topology

Buffer arrays for the left sided advective fluxes added in arrays_3d.
+flux_s_u, +flux_s_v, +flux_s_w, +diss_s_u, +diss_s_v, +diss_s_w,
+flux_s_pt, +diss_s_pt, +flux_s_e, +diss_s_e, +flux_s_q, +diss_s_q,
+flux_s_sa, +diss_s_sa
3d arrays for dissipation control added. (only necessary for vector arch.)
+var_x, +var_y, +var_z, +gamma_x, +gamma_y, +gamma_z
Default of momentum_advec and scalar_advec changed to 'ws-scheme' .
+exchange_mg added in control_parameters to steer the data exchange.
Parameters +nbgp, +nxlg, +nxrg, +nysg, +nyng added in indices.
flag array +boundary_flags added in indices to steer the degradation of order
of the advective fluxes when non-cyclic boundaries are used.
MPI-datatypes +type_y, +type_y_int and +type_yz for data_exchange added in
pegrid.
+sums_wsus_ws_l, +sums_wsvs_ws_l, +sums_us2_ws_l, +sums_vs2_ws_l,
+sums_ws2_ws_l, +sums_wspts_ws_l, +sums_wssas_ws_l, +sums_wsqs_ws_l
and +weight_substep added in statistics to steer the statistical evaluation
of turbulent fluxes in the advection routines.
LOGICALS +ws_scheme_sca and +ws_scheme_mom added to get a better performance
in prognostic_equations.
LOGICAL +dissipation_control control added to steer numerical dissipation
in ws-scheme.

Changed length of string run_description_header

pres.f90

New allocation of tend when ws-scheme and multigrid is used. This is due to
reasons of perforance of the data_exchange. The same is done with p after
poismg is called.
nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng when no
multigrid is used. Calls of exchange_horiz are modified.

bugfix: After pressure correction no volume flow correction in case of
non-cyclic boundary conditions
(has to be done only before pressure correction)

Call of SOR routine is referenced with ddzu_pres.

prognostic_equations.f90

Calls of the advection routines with WS5 added.
Calls of ws_statistics added to set the statistical arrays to zero after each
time step.

advec_particles.f90

Declaration of de_dx, de_dy, de_dz adapted to additional ghost points.
Furthermore the calls of exchange_horiz were modified.

asselin_filter.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

average_3d_data.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

boundary_conds.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
Removed mirror boundary conditions for u and v at the bottom in case of
ibc_uv_b == 0. Instead, dirichelt boundary conditions (u=v=0) are set
in init_3d_model

calc_liquid_water_content.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

calc_spectra.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng for
allocation of tend.

check_open.f90

Output of total array size was adapted to nbgp.

data_output_2d.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays local_2d and total_2d.
Calls of exchange_horiz are modified.

data_output_2d.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays. Calls of exchange_horiz are modified.
Skip-value skip_do_avs changed to a dynamic adaption of ghost points.

data_output_mask.f90

Calls of exchange_horiz are modified.

diffusion_e.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_s.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_u.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_v.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusion_w.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusivities.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

diffusivities.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Calls of exchange_horiz are modified.

exchange_horiz_2d.f90

Dynamic exchange of ghost points with nbgp, which depends on the advection
scheme. Exchange between left and right PEs is now done with MPI-vectors.

global_min_max.f90

Adapting of the index arrays, because MINLOC assumes lowerbound
at 1 and not at nbgp.

init_3d_model.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in loops and
allocation of arrays. Calls of exchange_horiz are modified.
Call ws_init to initialize arrays needed for statistical evaluation and
optimization when ws-scheme is used.
Initial volume flow is now calculated by using the variable hom_sum.
Therefore the correction of initial volume flow for non-flat topography
removed (removed u_nzb_p1_for_vfc and v_nzb_p1_for_vfc)
Changed surface boundary conditions for u and v in case of ibc_uv_b == 0 from
mirror bc to dirichlet boundary conditions (u=v=0), so that k=nzb is
representative for the height z0

Bugfix: type conversion of '1' to 64bit for the MAX function (ngp_3d_inner)

init_coupling.f90

determination of target_id's moved to init_pegrid

init_pt_anomaly.f90

Call of exchange_horiz are modified.

init_rankine.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Calls of exchange_horiz are modified.

init_slope.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

header.f90

Output of advection scheme.

poismg.f90

Calls of exchange_horiz are modified.

prandtl_fluxes.f90

Changed surface boundary conditions for u and v from mirror bc to dirichelt bc,
therefore u(uzb,:,:) and v(nzb,:,:) is now representative for the height z0
nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

production_e.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng

read_3d_binary.f90

+/- 1 replaced with +/- nbgp when swapping and allocating variables.

sor.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.
Call of exchange_horiz are modified.
bug removed in declaration of ddzw(), nz replaced by nzt+1

subsidence.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

sum_up_3d_data.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

surface_coupler.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng in
MPI_SEND() and MPI_RECV.
additional case for nonequivalent processor and grid topopolgy in ocean and
atmosphere added (coupling_topology = 1)

Added exchange of u and v from Ocean to Atmosphere

time_integration.f90

Calls of exchange_horiz are modified.
Adaption to slooping surface.

timestep.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

user_3d_data_averaging.f90, user_data_output_2d.f90, user_data_output_3d.f90,
user_actions.f90, user_init.f90, user_init_plant_canopy.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

user_read_restart_data.f90

Allocation with nbgp.

wall_fluxes.f90

nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng.

write_compressed.f90

Array bounds and nx, ny adapted with nbgp.

sor.f90

bug removed in declaration of ddzw(), nz replaced by nzt+1

Property svn:keywords set to Id

File size: 51.5 KB

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

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

Context Navigation

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

Download in other formats: