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

Last change on this file since 1007 was 1007, checked in by franke, 9 years ago

Bugfixes


Missing calculation of mean particle weighting factor for output added. (data_output_2d, data_output_3d, data_output_mask, sum_up_3d_data)
Calculation of mean particle radius for output now considers the weighting factor. (data_output_mask)
Calculation of sugrid-scale buoyancy flux for humidity and cloud droplets corrected. (flow_statistics)
Factor in calculation of enhancement factor for collision efficencies corrected. (lpm_collision_kernels)
Calculation of buoyancy production now considers the liquid water mixing ratio in case of cloud droplets. (production_e)

Changes


Calculation of buoyancy flux for humidity in case of WS-scheme is now using turbulent fluxes of WS-scheme. (flow_statistics)
Calculation of the collision kernels now in SI units. (lpm_collision_kernels)

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