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

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

optimization of two-moments cloud physics

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