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

Last change on this file since 667 was 667, checked in by suehring, 11 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 w*p* 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)       ! w*u*
161          sums_l(:,15,0) = sums_wsvs_ws_l(:,sr)       ! w*v*
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)      ! w*pts* from advec_s_ws
176          IF ( ocean ) sums_l(:,66,0) = sums_wssas_ws_l(:,sr) ! w*sa*
177          IF ( humidity  .OR.  passive_scalar ) sums_l(:,49,0) =              &
178                                                   sums_wsqs_ws_l(:,sr) !w*q*
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_cp*pt_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 w*e*
721!--             has to be adjusted
722                sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5 *          &
723                                             ( ust**2 + vst**2 + 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 w*u*
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 w*v*
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 u*pt*, v*pt*
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)     ! w*u*
991       hom(:,1,14,sr) = sums(:,14)     ! w"v"
992       hom(:,1,15,sr) = sums(:,15)     ! w*v*
993       hom(:,1,16,sr) = sums(:,16)     ! w"pt"
994       hom(:,1,17,sr) = sums(:,17)     ! w*pt*
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)     ! w*pt*BC
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)     ! w*2pt*
1009       hom(:,1,36,sr) = sums(:,36)     ! w*pt*2
1010       hom(:,1,37,sr) = sums(:,37)     ! w*e*
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)     ! w*vpt*       
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)     ! w*q* (w*qv*)
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)     ! w*qv*       
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)     ! w*u*u*/dz
1025       hom(:,1,56,sr) = sums(:,56)     ! w*p*/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)     ! u*pt*
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)     ! v*pt*
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)     ! w*sa*
1036       hom(:,1,67,sr) = sums(:,65) + sums(:,66)    ! wsa
1037       hom(:,1,68,sr) = sums(:,68)     ! w*p*
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.