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

Last change on this file since 1037 was 1037, checked in by raasch, 9 years ago

last commit documented

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