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

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

code has been put under the GNU General Public License (v3)

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