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

Last change on this file since 1116 was 1116, checked in by hoffmann, 8 years ago

last commit documented

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