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

Last change on this file since 1386 was 1386, checked in by boeske, 7 years ago

bugfix: simulated time before last timestep must be used for the profile data output of large scale forcing tendencies

  • Property svn:keywords set to Id
File size: 130.9 KB
Line 
1#if ! defined( __openacc )
2 SUBROUTINE flow_statistics
3
4!--------------------------------------------------------------------------------!
5! This file is part of PALM.
6!
7! PALM is free software: you can redistribute it and/or modify it under the terms
8! of the GNU General Public License as published by the Free Software Foundation,
9! either version 3 of the License, or (at your option) any later version.
10!
11! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
12! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
13! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
14!
15! You should have received a copy of the GNU General Public License along with
16! PALM. If not, see <http://www.gnu.org/licenses/>.
17!
18! Copyright 1997-2014 Leibniz Universitaet Hannover
19!--------------------------------------------------------------------------------!
20!
21! Current revisions:
22! -----------------
23! bugfix: simulated time before the last timestep is needed for the correct
24! calculation of the profiles of large scale forcing tendencies
25!
26! Former revisions:
27! -----------------
28! $Id: flow_statistics.f90 1386 2014-05-05 13:55:30Z boeske $
29!
30! 1382 2014-04-30 12:15:41Z boeske
31! Renamed variables which store large scale forcing tendencies
32! pt_lsa -> td_lsa_lpt, pt_subs -> td_sub_lpt,
33! q_lsa  -> td_lsa_q,   q_subs  -> td_sub_q,
34! added Neumann boundary conditions for profile data output of large scale
35! advection and subsidence terms at nzt+1
36!
37! 1374 2014-04-25 12:55:07Z raasch
38! bugfix: syntax errors removed from openacc-branch
39! missing variables added to ONLY-lists
40!
41! 1365 2014-04-22 15:03:56Z boeske
42! Output of large scale advection, large scale subsidence and nudging tendencies
43! +sums_ls_l, ngp_sums_ls, use_subsidence_tendencies
44!
45! 1353 2014-04-08 15:21:23Z heinze
46! REAL constants provided with KIND-attribute
47!
48! 1322 2014-03-20 16:38:49Z raasch
49! REAL constants defined as wp-kind
50!
51! 1320 2014-03-20 08:40:49Z raasch
52! ONLY-attribute added to USE-statements,
53! kind-parameters added to all INTEGER and REAL declaration statements,
54! kinds are defined in new module kinds,
55! revision history before 2012 removed,
56! comment fields (!:) to be used for variable explanations added to
57! all variable declaration statements
58!
59! 1299 2014-03-06 13:15:21Z heinze
60! Output of large scale vertical velocity w_subs
61!
62! 1257 2013-11-08 15:18:40Z raasch
63! openacc "end parallel" replaced by "end parallel loop"
64!
65! 1241 2013-10-30 11:36:58Z heinze
66! Output of ug and vg
67!
68! 1221 2013-09-10 08:59:13Z raasch
69! ported for openACC in separate #else branch
70!
71! 1179 2013-06-14 05:57:58Z raasch
72! comment for profile 77 added
73!
74! 1115 2013-03-26 18:16:16Z hoffmann
75! ql is calculated by calc_liquid_water_content
76!
77! 1111 2013-03-08 23:54:10Z raasch
78! openACC directive added
79!
80! 1053 2012-11-13 17:11:03Z hoffmann
81! additions for two-moment cloud physics scheme:
82! +nr, qr, qc, prr
83!
84! 1036 2012-10-22 13:43:42Z raasch
85! code put under GPL (PALM 3.9)
86!
87! 1007 2012-09-19 14:30:36Z franke
88! Calculation of buoyancy flux for humidity in case of WS-scheme is now using
89! turbulent fluxes of WS-scheme
90! Bugfix: Calculation of subgridscale buoyancy flux for humidity and cloud
91! droplets at nzb and nzt added
92!
93! 801 2012-01-10 17:30:36Z suehring
94! Calculation of turbulent fluxes in advec_ws is now thread-safe.
95!
96! Revision 1.1  1997/08/11 06:15:17  raasch
97! Initial revision
98!
99!
100! Description:
101! ------------
102! Compute average profiles and further average flow quantities for the different
103! user-defined (sub-)regions. The region indexed 0 is the total model domain.
104!
105! NOTE: For simplicity, nzb_s_inner and nzb_diff_s_inner are being used as a
106! ----  lower vertical index for k-loops for all variables, although strictly
107! speaking the k-loops would have to be split up according to the staggered
108! grid. However, this implies no error since staggered velocity components are
109! zero at the walls and inside buildings.
110!------------------------------------------------------------------------------!
111
112    USE arrays_3d,                                                             &
113        ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, p, prho, pt, q, qc, ql, qr, qs, &
114               qsws, qswst, rho, sa, saswsb, saswst, shf, td_lsa_lpt, td_lsa_q,&
115               td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us, usws,    &
116               uswst, vsws, v, vg, vpt, vswst, w, w_subs, zw
117       
118    USE cloud_parameters,                                                      &
119        ONLY :  l_d_cp, prr, pt_d_t
120       
121    USE control_parameters,                                                    &
122        ONLY :  average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
123                dt_3d, g, humidity, icloud_scheme, kappa, large_scale_forcing, &
124                large_scale_subsidence, max_pr_user, message_string, ocean,    &
125                passive_scalar, precipitation, simulated_time,                 &
126                use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
127                ws_scheme_mom, ws_scheme_sca
128       
129    USE cpulog,                                                                &
130        ONLY :  cpu_log, log_point
131       
132    USE grid_variables,                                                        &
133        ONLY :  ddx, ddy
134       
135    USE indices,                                                               &
136        ONLY :  ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums,      &
137                ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzb_diff_s_inner,        &
138                nzb_s_inner, nzt, nzt_diff
139       
140    USE kinds
141   
142    USE pegrid
143   
144    USE statistics
145
146    IMPLICIT NONE
147
148    INTEGER(iwp) ::  i                   !:
149    INTEGER(iwp) ::  j                   !:
150    INTEGER(iwp) ::  k                   !:
151    INTEGER(iwp) ::  nt                  !:
152    INTEGER(iwp) ::  omp_get_thread_num  !:
153    INTEGER(iwp) ::  sr                  !:
154    INTEGER(iwp) ::  tn                  !:
155   
156    LOGICAL ::  first  !:
157   
158    REAL(wp) ::  dptdz_threshold  !:
159    REAL(wp) ::  fac              !:
160    REAL(wp) ::  height           !:
161    REAL(wp) ::  pts              !:
162    REAL(wp) ::  sums_l_eper      !:
163    REAL(wp) ::  sums_l_etot      !:
164    REAL(wp) ::  ust              !:
165    REAL(wp) ::  ust2             !:
166    REAL(wp) ::  u2               !:
167    REAL(wp) ::  vst              !:
168    REAL(wp) ::  vst2             !:
169    REAL(wp) ::  v2               !:
170    REAL(wp) ::  w2               !:
171    REAL(wp) ::  z_i(2)           !:
172   
173    REAL(wp) ::  dptdz(nzb+1:nzt+1)    !:
174    REAL(wp) ::  sums_ll(nzb:nzt+1,2)  !:
175
176    CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
177
178    !$acc update host( km, kh, e, pt, qs, qsws, rif, shf, ts, u, usws, v, vsws, w )
179
180!
181!-- To be on the safe side, check whether flow_statistics has already been
182!-- called once after the current time step
183    IF ( flow_statistics_called )  THEN
184
185       message_string = 'flow_statistics is called two times within one ' // &
186                        'timestep'
187       CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
188
189    ENDIF
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_wp
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_wp *                        & 
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_wp * ( &
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_wp
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_wp * ( 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_wp
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_wp *                &
545                                  ( ust2 + vst2 + w2 ) * rmask(j,i,sr)
546                   sums_l_eper     = sums_l_eper +                             &
547                                     0.5_wp * ( 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!
557!--    Horizontally averaged profiles of the vertical fluxes
558
559       !$OMP DO
560       DO  i = nxl, nxr
561          DO  j = nys, nyn
562!
563!--          Subgridscale fluxes (without Prandtl layer from k=nzb,
564!--          oterwise from k=nzb+1)
565!--          NOTE: for simplicity, nzb_diff_s_inner is used below, although
566!--          ----  strictly speaking the following k-loop would have to be
567!--                split up according to the staggered grid.
568!--                However, this implies no error since staggered velocity
569!--                components are zero at the walls and inside buildings.
570
571             DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
572!
573!--             Momentum flux w"u"
574                sums_l(k,12,tn) = sums_l(k,12,tn) - 0.25_wp * (                &
575                               km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
576                                                           ) * (               &
577                                   ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
578                                 + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
579                                                               ) * rmask(j,i,sr)
580!
581!--             Momentum flux w"v"
582                sums_l(k,14,tn) = sums_l(k,14,tn) - 0.25_wp * (                &
583                               km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
584                                                           ) * (               &
585                                   ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
586                                 + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
587                                                               ) * rmask(j,i,sr)
588!
589!--             Heat flux w"pt"
590                sums_l(k,16,tn) = sums_l(k,16,tn)                              &
591                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
592                                               * ( pt(k+1,j,i) - pt(k,j,i) )   &
593                                               * ddzu(k+1) * rmask(j,i,sr)
594
595
596!
597!--             Salinity flux w"sa"
598                IF ( ocean )  THEN
599                   sums_l(k,65,tn) = sums_l(k,65,tn)                           &
600                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
601                                               * ( sa(k+1,j,i) - sa(k,j,i) )   &
602                                               * ddzu(k+1) * rmask(j,i,sr)
603                ENDIF
604
605!
606!--             Buoyancy flux, water flux (humidity flux) w"q"
607                IF ( humidity ) THEN
608                   sums_l(k,45,tn) = sums_l(k,45,tn)                           &
609                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
610                                               * ( vpt(k+1,j,i) - vpt(k,j,i) ) &
611                                               * ddzu(k+1) * rmask(j,i,sr)
612                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
613                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
614                                               * ( q(k+1,j,i) - q(k,j,i) )     &
615                                               * ddzu(k+1) * rmask(j,i,sr)
616
617                   IF ( cloud_physics ) THEN
618                      sums_l(k,51,tn) = sums_l(k,51,tn)                        &
619                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
620                                               * ( ( q(k+1,j,i) - ql(k+1,j,i) )&
621                                                - ( q(k,j,i) - ql(k,j,i) ) )   &
622                                               * ddzu(k+1) * rmask(j,i,sr) 
623                   ENDIF
624                ENDIF
625
626!
627!--             Passive scalar flux
628                IF ( passive_scalar )  THEN
629                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
630                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
631                                               * ( q(k+1,j,i) - q(k,j,i) )     &
632                                               * ddzu(k+1) * rmask(j,i,sr)
633                ENDIF
634
635             ENDDO
636
637!
638!--          Subgridscale fluxes in the Prandtl layer
639             IF ( use_surface_fluxes )  THEN
640                sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
641                                    usws(j,i) * rmask(j,i,sr)     ! w"u"
642                sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
643                                    vsws(j,i) * rmask(j,i,sr)     ! w"v"
644                sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
645                                    shf(j,i)  * rmask(j,i,sr)     ! w"pt"
646                sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
647                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
648                sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
649                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
650                IF ( ocean )  THEN
651                   sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
652                                       saswsb(j,i) * rmask(j,i,sr)  ! w"sa"
653                ENDIF
654                IF ( humidity )  THEN
655                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
656                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
657                   sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                   &
658                                       ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) *     &
659                                       shf(j,i) + 0.61_wp * pt(nzb,j,i) *      &
660                                                  qsws(j,i) )
661                   IF ( cloud_droplets )  THEN
662                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                &
663                                         ( 1.0_wp + 0.61_wp * q(nzb,j,i) -     &
664                                           ql(nzb,j,i) ) * shf(j,i) +          &
665                                           0.61_wp * pt(nzb,j,i) * qsws(j,i) )
666                   ENDIF
667                   IF ( cloud_physics )  THEN
668!
669!--                   Formula does not work if ql(nzb) /= 0.0
670                      sums_l(nzb,51,tn) = sums_l(nzb,51,tn) + &   ! w"q" (w"qv")
671                                          qsws(j,i) * rmask(j,i,sr)
672                   ENDIF
673                ENDIF
674                IF ( passive_scalar )  THEN
675                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
676                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
677                ENDIF
678             ENDIF
679
680!
681!--          Subgridscale fluxes at the top surface
682             IF ( use_top_fluxes )  THEN
683                sums_l(nzt:nzt+1,12,tn) = sums_l(nzt:nzt+1,12,tn) + &
684                                    uswst(j,i) * rmask(j,i,sr)    ! w"u"
685                sums_l(nzt:nzt+1,14,tn) = sums_l(nzt:nzt+1,14,tn) + &
686                                    vswst(j,i) * rmask(j,i,sr)    ! w"v"
687                sums_l(nzt:nzt+1,16,tn) = sums_l(nzt:nzt+1,16,tn) + &
688                                    tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
689                sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + &
690                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
691                sums_l(nzt:nzt+1,61,tn) = sums_l(nzt:nzt+1,61,tn) + &
692                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
693
694                IF ( ocean )  THEN
695                   sums_l(nzt,65,tn) = sums_l(nzt,65,tn) + &
696                                       saswst(j,i) * rmask(j,i,sr)  ! w"sa"
697                ENDIF
698                IF ( humidity )  THEN
699                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) +                     &
700                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
701                   sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                   &
702                                       ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) *     &
703                                       tswst(j,i) + 0.61_wp * pt(nzt,j,i) *    &
704                                                             qswst(j,i) )
705                   IF ( cloud_droplets )  THEN
706                      sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                &
707                                          ( 1.0_wp + 0.61_wp * q(nzt,j,i) -    &
708                                            ql(nzt,j,i) ) * tswst(j,i) +       &
709                                            0.61_wp * pt(nzt,j,i) * qswst(j,i) )
710                   ENDIF
711                   IF ( cloud_physics )  THEN
712!
713!--                   Formula does not work if ql(nzb) /= 0.0
714                      sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + &   ! w"q" (w"qv")
715                                          qswst(j,i) * rmask(j,i,sr)
716                   ENDIF
717                ENDIF
718                IF ( passive_scalar )  THEN
719                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
720                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
721                ENDIF
722             ENDIF
723
724!
725!--          Resolved fluxes (can be computed for all horizontal points)
726!--          NOTE: for simplicity, nzb_s_inner is used below, although strictly
727!--          ----  speaking the following k-loop would have to be split up and
728!--                rearranged according to the staggered grid.
729             DO  k = nzb_s_inner(j,i), nzt
730                ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +                  &
731                                 u(k+1,j,i) - hom(k+1,1,1,sr) )
732                vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +                  &
733                                 v(k+1,j,i) - hom(k+1,1,2,sr) )
734                pts = 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) +                 &
735                                 pt(k+1,j,i) - hom(k+1,1,4,sr) )
736
737!--             Higher moments
738                sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 *        &
739                                                    rmask(j,i,sr)
740                sums_l(k,36,tn) = sums_l(k,36,tn) + pts**2 * w(k,j,i) *        &
741                                                    rmask(j,i,sr)
742
743!
744!--             Salinity flux and density (density does not belong to here,
745!--             but so far there is no other suitable place to calculate)
746                IF ( ocean )  THEN
747                   IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
748                      pts = 0.5_wp * ( sa(k,j,i)   - hom(k,1,23,sr) +          &
749                                       sa(k+1,j,i) - hom(k+1,1,23,sr) )
750                      sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) *     &
751                                        rmask(j,i,sr)
752                   ENDIF
753                   sums_l(k,64,tn) = sums_l(k,64,tn) + rho(k,j,i) *            &
754                                                       rmask(j,i,sr)
755                   sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) *           &
756                                                       rmask(j,i,sr)
757                ENDIF
758
759!
760!--             Buoyancy flux, water flux, humidity flux, liquid water
761!--             content, rain drop concentration and rain water content
762                IF ( humidity )  THEN
763                   IF ( cloud_physics .OR. cloud_droplets )  THEN
764                      pts = 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +         &
765                                    vpt(k+1,j,i) - hom(k+1,1,44,sr) )
766                      sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *     &
767                                                          rmask(j,i,sr)
768                      IF ( .NOT. cloud_droplets )  THEN
769                         pts = 0.5_wp *                                        &
770                              ( ( q(k,j,i) - ql(k,j,i) ) -                     &
771                              hom(k,1,42,sr) +                                 &
772                              ( q(k+1,j,i) - ql(k+1,j,i) ) -                   &
773                              hom(k+1,1,42,sr) )
774                         sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) *  &
775                                                             rmask(j,i,sr)
776                         IF ( icloud_scheme == 0  )  THEN
777                            sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) *    &
778                                                                rmask(j,i,sr)
779                            sums_l(k,75,tn) = sums_l(k,75,tn) + qc(k,j,i) *    &
780                                                                rmask(j,i,sr)
781                            IF ( precipitation )  THEN
782                               sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) * &
783                                                                   rmask(j,i,sr)
784                               sums_l(k,74,tn) = sums_l(k,74,tn) + qr(k,j,i) * &
785                                                                   rmask(j,i,sr)
786                               sums_l(k,76,tn) = sums_l(k,76,tn) + prr(k,j,i) *&
787                                                                   rmask(j,i,sr)
788                            ENDIF
789                         ELSE
790                            sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) *    &
791                                                                rmask(j,i,sr)
792                         ENDIF
793                      ELSE
794                         sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) *       &
795                                                             rmask(j,i,sr)
796                      ENDIF
797                   ELSE
798                      IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
799                         pts = 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +      &
800                                          vpt(k+1,j,i) - hom(k+1,1,44,sr) )
801                         sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *  &
802                                                             rmask(j,i,sr)
803                      ELSE IF ( ws_scheme_sca .AND. sr == 0 )  THEN
804                         sums_l(k,46,tn) = ( 1.0_wp + 0.61_wp *                & 
805                                             hom(k,1,41,sr) ) *                &
806                                           sums_l(k,17,tn) +                   &
807                                           0.61_wp * hom(k,1,4,sr) *           &
808                                           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_wp * ( 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_wp *        &
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_wp * ( 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_wp * ( 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_wp *                 &
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_wp *                 &
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_wp *                 &
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_wp * ( 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!--    Density at top follows Neumann condition
884       IF ( ocean )  THEN
885          sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
886          sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
887       ENDIF
888
889!
890!--    Divergence of vertical flux of resolved scale energy and pressure
891!--    fluctuations as well as flux of pressure fluctuation itself (68).
892!--    First calculate the products, then the divergence.
893!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
894       IF ( hom(nzb+1,2,55,0) /= 0.0_wp  .OR.  hom(nzb+1,2,68,0) /= 0.0_wp )  THEN
895
896          sums_ll = 0.0_wp  ! local array
897
898          !$OMP DO
899          DO  i = nxl, nxr
900             DO  j = nys, nyn
901                DO  k = nzb_s_inner(j,i)+1, nzt
902
903                   sums_ll(k,1) = sums_ll(k,1) + 0.5_wp * w(k,j,i) * (         &
904                  ( 0.25_wp * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1)    &
905                            - 0.5_wp * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) )     &
906                              ) )**2                                           &
907                + ( 0.25_wp * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i)    &
908                            - 0.5_wp * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) )     &
909                              ) )**2                                           &
910                + w(k,j,i)**2                                        )
911
912                   sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i)             &
913                                               * ( p(k,j,i) + p(k+1,j,i) )
914
915                ENDDO
916             ENDDO
917          ENDDO
918          sums_ll(0,1)     = 0.0_wp    ! because w is zero at the bottom
919          sums_ll(nzt+1,1) = 0.0_wp
920          sums_ll(0,2)     = 0.0_wp
921          sums_ll(nzt+1,2) = 0.0_wp
922
923          DO  k = nzb+1, nzt
924             sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
925             sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
926             sums_l(k,68,tn) = sums_ll(k,2)
927          ENDDO
928          sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
929          sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
930          sums_l(nzb,68,tn) = 0.0_wp    ! because w* = 0 at nzb
931
932       ENDIF
933
934!
935!--    Divergence of vertical flux of SGS TKE and the flux itself (69)
936       IF ( hom(nzb+1,2,57,0) /= 0.0_wp  .OR.  hom(nzb+1,2,69,0) /= 0.0_wp )  THEN
937
938          !$OMP DO
939          DO  i = nxl, nxr
940             DO  j = nys, nyn
941                DO  k = nzb_s_inner(j,i)+1, nzt
942
943                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5_wp * (              &
944                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
945                 - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
946                                                                ) * ddzw(k)
947
948                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5_wp * (              &
949                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
950                                                                )
951
952                ENDDO
953             ENDDO
954          ENDDO
955          sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
956          sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
957
958       ENDIF
959
960!
961!--    Horizontal heat fluxes (subgrid, resolved, total).
962!--    Do it only, if profiles shall be plotted.
963       IF ( hom(nzb+1,2,58,0) /= 0.0_wp ) THEN
964
965          !$OMP DO
966          DO  i = nxl, nxr
967             DO  j = nys, nyn
968                DO  k = nzb_s_inner(j,i)+1, nzt
969!
970!--                Subgrid horizontal heat fluxes u"pt", v"pt"
971                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5_wp *                &
972                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
973                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
974                                                 * ddx * rmask(j,i,sr)
975                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
976                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
977                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
978                                                 * ddy * rmask(j,i,sr)
979!
980!--                Resolved horizontal heat fluxes u*pt*, v*pt*
981                   sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
982                                                  ( u(k,j,i) - hom(k,1,1,sr) ) &
983                                    * 0.5_wp * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
984                                                 pt(k,j,i)   - hom(k,1,4,sr) )
985                   pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
986                                    pt(k,j,i)   - hom(k,1,4,sr) )
987                   sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
988                                                  ( v(k,j,i) - hom(k,1,2,sr) ) &
989                                    * 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
990                                                 pt(k,j,i)   - hom(k,1,4,sr) )
991                ENDDO
992             ENDDO
993          ENDDO
994!
995!--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
996          sums_l(nzb,58,tn) = 0.0_wp
997          sums_l(nzb,59,tn) = 0.0_wp
998          sums_l(nzb,60,tn) = 0.0_wp
999          sums_l(nzb,61,tn) = 0.0_wp
1000          sums_l(nzb,62,tn) = 0.0_wp
1001          sums_l(nzb,63,tn) = 0.0_wp
1002
1003       ENDIF
1004
1005!
1006!--    Collect current large scale advection and subsidence tendencies for
1007!--    data output
1008       IF ( large_scale_forcing )  THEN
1009!
1010!--       Interpolation in time of LSF_DATA
1011          nt = 1
1012          DO WHILE ( simulated_time - dt_3d > time_vert(nt) )
1013             nt = nt + 1
1014          ENDDO
1015          IF ( simulated_time - dt_3d /= time_vert(nt) )  THEN
1016            nt = nt - 1
1017          ENDIF
1018
1019          fac = ( simulated_time - dt_3d - time_vert(nt) )                     &
1020                / ( time_vert(nt+1)-time_vert(nt) )
1021
1022
1023          DO  k = nzb, nzt
1024             sums_ls_l(k,0) = td_lsa_lpt(k,nt)                                 &
1025                              + fac * ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )
1026             sums_ls_l(k,1) = td_lsa_q(k,nt)                                   &
1027                              + fac * ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) )
1028          ENDDO
1029
1030          sums_ls_l(nzt+1,0) = sums_ls_l(nzt,0)
1031          sums_ls_l(nzt+1,1) = sums_ls_l(nzt,1)
1032
1033          IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
1034
1035             DO  k = nzb, nzt
1036                sums_ls_l(k,2) = td_sub_lpt(k,nt) + fac *                      &
1037                                 ( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )
1038                sums_ls_l(k,3) = td_sub_q(k,nt) + fac *                        &
1039                                 ( td_sub_q(k,nt+1) - td_sub_q(k,nt) )
1040             ENDDO
1041
1042             sums_ls_l(nzt+1,2) = sums_ls_l(nzt,2)
1043             sums_ls_l(nzt+1,3) = sums_ls_l(nzt,3)
1044
1045          ENDIF
1046
1047       ENDIF
1048
1049!
1050!--    Calculate the user-defined profiles
1051       CALL user_statistics( 'profiles', sr, tn )
1052       !$OMP END PARALLEL
1053
1054!
1055!--    Summation of thread sums
1056       IF ( threads_per_task > 1 )  THEN
1057          DO  i = 1, threads_per_task-1
1058             sums_l(:,3,0)          = sums_l(:,3,0) + sums_l(:,3,i)
1059             sums_l(:,4:40,0)       = sums_l(:,4:40,0) + sums_l(:,4:40,i)
1060             sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
1061                                      sums_l(:,45:pr_palm,i)
1062             IF ( max_pr_user > 0 )  THEN
1063                sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
1064                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
1065                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
1066             ENDIF
1067          ENDDO
1068       ENDIF
1069
1070#if defined( __parallel )
1071
1072!
1073!--    Compute total sum from local sums
1074       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1075       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL,   &
1076                           MPI_SUM, comm2d, ierr )
1077       IF ( large_scale_forcing )  THEN
1078          CALL MPI_ALLREDUCE( sums_ls_l(nzb,2), sums(nzb,83), ngp_sums_ls,     &
1079                              MPI_REAL, MPI_SUM, comm2d, ierr )
1080       ENDIF
1081#else
1082       sums = sums_l(:,:,0)
1083       IF ( large_scale_forcing )  THEN
1084          sums(:,81:88) = sums_ls_l
1085       ENDIF
1086#endif
1087
1088!
1089!--    Final values are obtained by division by the total number of grid points
1090!--    used for summation. After that store profiles.
1091!--    Profiles:
1092       DO  k = nzb, nzt+1
1093          sums(k,3)               = sums(k,3)           / ngp_2dh(sr)
1094          sums(k,8:11)            = sums(k,8:11)        / ngp_2dh_s_inner(k,sr)
1095          sums(k,12:22)           = sums(k,12:22)       / ngp_2dh(sr)
1096          sums(k,23:29)           = sums(k,23:29)       / ngp_2dh_s_inner(k,sr)
1097          sums(k,30:32)           = sums(k,30:32)       / ngp_2dh(sr)
1098          sums(k,33:34)           = sums(k,33:34)       / ngp_2dh_s_inner(k,sr)
1099          sums(k,35:39)           = sums(k,35:39)       / ngp_2dh(sr)
1100          sums(k,40)              = sums(k,40)          / ngp_2dh_s_inner(k,sr)
1101          sums(k,45:53)           = sums(k,45:53)       / ngp_2dh(sr)
1102          sums(k,54)              = sums(k,54)          / ngp_2dh_s_inner(k,sr)
1103          sums(k,55:63)           = sums(k,55:63)       / ngp_2dh(sr)
1104          sums(k,64)              = sums(k,64)          / ngp_2dh_s_inner(k,sr)
1105          sums(k,65:69)           = sums(k,65:69)       / ngp_2dh(sr)
1106          sums(k,70:80)           = sums(k,70:80)       / ngp_2dh_s_inner(k,sr)
1107          sums(k,81:88)           = sums(k,81:88)       / ngp_2dh(sr)
1108          sums(k,89:pr_palm-2)    = sums(k,89:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
1109       ENDDO
1110
1111!--    Upstream-parts
1112       sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
1113!--    u* and so on
1114!--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
1115!--    size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
1116!--    above the topography, they are being divided by ngp_2dh(sr)
1117       sums(nzb:nzb+3,pr_palm)    = sums(nzb:nzb+3,pr_palm)    / &
1118                                    ngp_2dh(sr)
1119       sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)       / &    ! qs
1120                                    ngp_2dh(sr)
1121!--    eges, e*
1122       sums(nzb+4:nzb+5,pr_palm)  = sums(nzb+4:nzb+5,pr_palm)  / &
1123                                    ngp_3d(sr)
1124!--    Old and new divergence
1125       sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
1126                                    ngp_3d_inner(sr)
1127
1128!--    User-defined profiles
1129       IF ( max_pr_user > 0 )  THEN
1130          DO  k = nzb, nzt+1
1131             sums(k,pr_palm+1:pr_palm+max_pr_user) = &
1132                                    sums(k,pr_palm+1:pr_palm+max_pr_user) / &
1133                                    ngp_2dh_s_inner(k,sr)
1134          ENDDO
1135       ENDIF
1136
1137!
1138!--    Collect horizontal average in hom.
1139!--    Compute deduced averages (e.g. total heat flux)
1140       hom(:,1,3,sr)  = sums(:,3)      ! w
1141       hom(:,1,8,sr)  = sums(:,8)      ! e     profiles 5-7 are initial profiles
1142       hom(:,1,9,sr)  = sums(:,9)      ! km
1143       hom(:,1,10,sr) = sums(:,10)     ! kh
1144       hom(:,1,11,sr) = sums(:,11)     ! l
1145       hom(:,1,12,sr) = sums(:,12)     ! w"u"
1146       hom(:,1,13,sr) = sums(:,13)     ! w*u*
1147       hom(:,1,14,sr) = sums(:,14)     ! w"v"
1148       hom(:,1,15,sr) = sums(:,15)     ! w*v*
1149       hom(:,1,16,sr) = sums(:,16)     ! w"pt"
1150       hom(:,1,17,sr) = sums(:,17)     ! w*pt*
1151       hom(:,1,18,sr) = sums(:,16) + sums(:,17)    ! wpt
1152       hom(:,1,19,sr) = sums(:,12) + sums(:,13)    ! wu
1153       hom(:,1,20,sr) = sums(:,14) + sums(:,15)    ! wv
1154       hom(:,1,21,sr) = sums(:,21)     ! w*pt*BC
1155       hom(:,1,22,sr) = sums(:,16) + sums(:,21)    ! wptBC
1156                                       ! profile 24 is initial profile (sa)
1157                                       ! profiles 25-29 left empty for initial
1158                                       ! profiles
1159       hom(:,1,30,sr) = sums(:,30)     ! u*2
1160       hom(:,1,31,sr) = sums(:,31)     ! v*2
1161       hom(:,1,32,sr) = sums(:,32)     ! w*2
1162       hom(:,1,33,sr) = sums(:,33)     ! pt*2
1163       hom(:,1,34,sr) = sums(:,34)     ! e*
1164       hom(:,1,35,sr) = sums(:,35)     ! w*2pt*
1165       hom(:,1,36,sr) = sums(:,36)     ! w*pt*2
1166       hom(:,1,37,sr) = sums(:,37)     ! w*e*
1167       hom(:,1,38,sr) = sums(:,38)     ! w*3
1168       hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20_wp )**1.5_wp   ! Sw
1169       hom(:,1,40,sr) = sums(:,40)     ! p
1170       hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
1171       hom(:,1,46,sr) = sums(:,46)     ! w*vpt*       
1172       hom(:,1,47,sr) = sums(:,45) + sums(:,46)    ! wvpt
1173       hom(:,1,48,sr) = sums(:,48)     ! w"q" (w"qv")
1174       hom(:,1,49,sr) = sums(:,49)     ! w*q* (w*qv*)
1175       hom(:,1,50,sr) = sums(:,48) + sums(:,49)    ! wq (wqv)
1176       hom(:,1,51,sr) = sums(:,51)     ! w"qv"
1177       hom(:,1,52,sr) = sums(:,52)     ! w*qv*       
1178       hom(:,1,53,sr) = sums(:,52) + sums(:,51)    ! wq (wqv)
1179       hom(:,1,54,sr) = sums(:,54)     ! ql
1180       hom(:,1,55,sr) = sums(:,55)     ! w*u*u*/dz
1181       hom(:,1,56,sr) = sums(:,56)     ! w*p*/dz
1182       hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho )/dz
1183       hom(:,1,58,sr) = sums(:,58)     ! u"pt"
1184       hom(:,1,59,sr) = sums(:,59)     ! u*pt*
1185       hom(:,1,60,sr) = sums(:,58) + sums(:,59)    ! upt_t
1186       hom(:,1,61,sr) = sums(:,61)     ! v"pt"
1187       hom(:,1,62,sr) = sums(:,62)     ! v*pt*
1188       hom(:,1,63,sr) = sums(:,61) + sums(:,62)    ! vpt_t
1189       hom(:,1,64,sr) = sums(:,64)     ! rho
1190       hom(:,1,65,sr) = sums(:,65)     ! w"sa"
1191       hom(:,1,66,sr) = sums(:,66)     ! w*sa*
1192       hom(:,1,67,sr) = sums(:,65) + sums(:,66)    ! wsa
1193       hom(:,1,68,sr) = sums(:,68)     ! w*p*
1194       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho
1195       hom(:,1,70,sr) = sums(:,70)     ! q*2
1196       hom(:,1,71,sr) = sums(:,71)     ! prho
1197       hom(:,1,72,sr) = hyp * 1E-4_wp  ! hyp in dbar
1198       hom(:,1,73,sr) = sums(:,73)     ! nr
1199       hom(:,1,74,sr) = sums(:,74)     ! qr
1200       hom(:,1,75,sr) = sums(:,75)     ! qc
1201       hom(:,1,76,sr) = sums(:,76)     ! prr (precipitation rate)
1202                                       ! 77 is initial density profile
1203       hom(:,1,78,sr) = ug             ! ug
1204       hom(:,1,79,sr) = vg             ! vg
1205       hom(:,1,80,sr) = w_subs         ! w_subs
1206
1207       IF ( large_scale_forcing )  THEN
1208          hom(:,1,81,sr) = sums_ls_l(:,0)          ! td_lsa_lpt
1209          hom(:,1,82,sr) = sums_ls_l(:,1)          ! td_lsa_q
1210          IF ( use_subsidence_tendencies )  THEN
1211             hom(:,1,83,sr) = sums_ls_l(:,2)       ! td_sub_lpt
1212             hom(:,1,84,sr) = sums_ls_l(:,3)       ! td_sub_q
1213          ELSE
1214             hom(:,1,83,sr) = sums(:,83)           ! td_sub_lpt
1215             hom(:,1,84,sr) = sums(:,84)           ! td_sub_q
1216          ENDIF
1217          hom(:,1,85,sr) = sums(:,85)              ! td_nud_lpt
1218          hom(:,1,86,sr) = sums(:,86)              ! td_nud_q
1219          hom(:,1,87,sr) = sums(:,87)              ! td_nud_u
1220          hom(:,1,88,sr) = sums(:,88)              ! td_nud_v
1221       ENDIF
1222
1223       hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
1224                                       ! upstream-parts u_x, u_y, u_z, v_x,
1225                                       ! v_y, usw. (in last but one profile)
1226       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
1227                                       ! u*, w'u', w'v', t* (in last profile)
1228
1229       IF ( max_pr_user > 0 )  THEN    ! user-defined profiles
1230          hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
1231                               sums(:,pr_palm+1:pr_palm+max_pr_user)
1232       ENDIF
1233
1234!
1235!--    Determine the boundary layer height using two different schemes.
1236!--    First scheme: Starting from the Earth's (Ocean's) surface, look for the
1237!--    first relative minimum (maximum) of the total heat flux.
1238!--    The corresponding height is assumed as the boundary layer height, if it
1239!--    is less than 1.5 times the height where the heat flux becomes negative
1240!--    (positive) for the first time.
1241       z_i(1) = 0.0_wp
1242       first = .TRUE.
1243
1244       IF ( ocean )  THEN
1245          DO  k = nzt, nzb+1, -1
1246             IF ( first .AND. hom(k,1,18,sr) < 0.0_wp                          &
1247                .AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp)  THEN
1248                first = .FALSE.
1249                height = zw(k)
1250             ENDIF
1251             IF ( hom(k,1,18,sr) < 0.0_wp  .AND.                               &
1252                  abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND.                        &
1253                  hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
1254                IF ( zw(k) < 1.5_wp * height )  THEN
1255                   z_i(1) = zw(k)
1256                ELSE
1257                   z_i(1) = height
1258                ENDIF
1259                EXIT
1260             ENDIF
1261          ENDDO
1262       ELSE
1263          DO  k = nzb, nzt-1
1264             IF ( first .AND. hom(k,1,18,sr) < 0.0_wp                          &
1265                .AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp )  THEN
1266                first = .FALSE.
1267                height = zw(k)
1268             ENDIF
1269             IF ( hom(k,1,18,sr) < 0.0_wp  .AND.                               & 
1270                  abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND.                        &
1271                  hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
1272                IF ( zw(k) < 1.5_wp * height )  THEN
1273                   z_i(1) = zw(k)
1274                ELSE
1275                   z_i(1) = height
1276                ENDIF
1277                EXIT
1278             ENDIF
1279          ENDDO
1280       ENDIF
1281
1282!
1283!--    Second scheme: Gradient scheme from Sullivan et al. (1998), modified
1284!--    by Uhlenbrock(2006). The boundary layer height is the height with the
1285!--    maximal local temperature gradient: starting from the second (the last
1286!--    but one) vertical gridpoint, the local gradient must be at least
1287!--    0.2K/100m and greater than the next four gradients.
1288!--    WARNING: The threshold value of 0.2K/100m must be adjusted for the
1289!--             ocean case!
1290       z_i(2) = 0.0_wp
1291       DO  k = nzb+1, nzt+1
1292          dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
1293       ENDDO
1294       dptdz_threshold = 0.2_wp / 100.0_wp
1295
1296       IF ( ocean )  THEN
1297          DO  k = nzt+1, nzb+5, -1
1298             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
1299                  dptdz(k) > dptdz(k-1)  .AND.  dptdz(k) > dptdz(k-2)  .AND.  &
1300                  dptdz(k) > dptdz(k-3)  .AND.  dptdz(k) > dptdz(k-4) )  THEN
1301                z_i(2) = zw(k-1)
1302                EXIT
1303             ENDIF
1304          ENDDO
1305       ELSE
1306          DO  k = nzb+1, nzt-3
1307             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
1308                  dptdz(k) > dptdz(k+1)  .AND.  dptdz(k) > dptdz(k+2)  .AND.  &
1309                  dptdz(k) > dptdz(k+3)  .AND.  dptdz(k) > dptdz(k+4) )  THEN
1310                z_i(2) = zw(k-1)
1311                EXIT
1312             ENDIF
1313          ENDDO
1314       ENDIF
1315
1316       hom(nzb+6,1,pr_palm,sr) = z_i(1)
1317       hom(nzb+7,1,pr_palm,sr) = z_i(2)
1318
1319!
1320!--    Computation of both the characteristic vertical velocity and
1321!--    the characteristic convective boundary layer temperature.
1322!--    The horizontal average at nzb+1 is input for the average temperature.
1323       IF ( hom(nzb,1,18,sr) > 0.0_wp .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8_wp  &
1324           .AND.  z_i(1) /= 0.0_wp )  THEN
1325          hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) *                 &
1326                                       hom(nzb,1,18,sr) *                      &
1327                                       ABS( z_i(1) ) )**0.333333333_wp
1328!--       so far this only works if Prandtl layer is used
1329          hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
1330       ELSE
1331          hom(nzb+8,1,pr_palm,sr)  = 0.0_wp
1332          hom(nzb+11,1,pr_palm,sr) = 0.0_wp
1333       ENDIF
1334
1335!
1336!--    Collect the time series quantities
1337       ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)     ! E
1338       ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)     ! E*
1339       ts_value(3,sr) = dt_3d
1340       ts_value(4,sr) = hom(nzb,1,pr_palm,sr)       ! u*
1341       ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr)     ! th*
1342       ts_value(6,sr) = u_max
1343       ts_value(7,sr) = v_max
1344       ts_value(8,sr) = w_max
1345       ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr)    ! new divergence
1346       ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr)    ! old Divergence
1347       ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr)    ! z_i(1)
1348       ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr)    ! z_i(2)
1349       ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr)    ! w*
1350       ts_value(14,sr) = hom(nzb,1,16,sr)           ! w'pt'   at k=0
1351       ts_value(15,sr) = hom(nzb+1,1,16,sr)         ! w'pt'   at k=1
1352       ts_value(16,sr) = hom(nzb+1,1,18,sr)         ! wpt     at k=1
1353       ts_value(17,sr) = hom(nzb,1,4,sr)            ! pt(0)
1354       ts_value(18,sr) = hom(nzb+1,1,4,sr)          ! pt(zp)
1355       ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr)    ! u'w'    at k=0
1356       ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)    ! v'w'    at k=0
1357       ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
1358
1359       IF ( ts_value(5,sr) /= 0.0_wp )  THEN
1360          ts_value(22,sr) = ts_value(4,sr)**2 / &
1361                            ( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
1362       ELSE
1363          ts_value(22,sr) = 10000.0_wp
1364       ENDIF
1365
1366       ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
1367!
1368!--    Calculate additional statistics provided by the user interface
1369       CALL user_statistics( 'time_series', sr, 0 )
1370
1371    ENDDO    ! loop of the subregions
1372
1373!
1374!-- If required, sum up horizontal averages for subsequent time averaging
1375    IF ( do_sum )  THEN
1376       IF ( average_count_pr == 0 )  hom_sum = 0.0_wp
1377       hom_sum = hom_sum + hom(:,1,:,:)
1378       average_count_pr = average_count_pr + 1
1379       do_sum = .FALSE.
1380    ENDIF
1381
1382!
1383!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
1384!-- This flag is reset after each time step in time_integration.
1385    flow_statistics_called = .TRUE.
1386
1387    CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
1388
1389
1390 END SUBROUTINE flow_statistics
1391
1392
1393#else
1394
1395
1396!------------------------------------------------------------------------------!
1397! flow statistics - accelerator version
1398!------------------------------------------------------------------------------!
1399 SUBROUTINE flow_statistics
1400
1401    USE arrays_3d,                                                             &
1402        ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, p, prho, pt, q, qc, ql, qr, qs, &
1403               qsws, qswst, rho, sa, saswsb, saswst, shf, td_lsa_lpt, td_lsa_q,&
1404               td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us, usws,    &
1405               uswst, vsws, v, vg, vpt, vswst, w, w_subs, zw
1406                 
1407       
1408    USE cloud_parameters,                                                      &
1409        ONLY:  l_d_cp, prr, pt_d_t
1410       
1411    USE control_parameters,                                                    &
1412        ONLY :  average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
1413                dt_3d, g, humidity, icloud_scheme, kappa, large_scale_forcing, &
1414                large_scale_subsidence, max_pr_user, message_string, ocean,    &
1415                passive_scalar, precipitation, simulated_time,                 &
1416                use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
1417                ws_scheme_mom, ws_scheme_sca
1418       
1419    USE cpulog,                                                                &
1420        ONLY:  cpu_log, log_point
1421       
1422    USE grid_variables,                                                        &
1423        ONLY:  ddx, ddy
1424       
1425    USE indices,                                                               &
1426        ONLY:  ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums, nxl,  &
1427               nxr, nyn, nys, nzb, nzb_diff_s_inner, nzb_s_inner, nzt,         &
1428               nzt_diff, rflags_invers
1429       
1430    USE kinds
1431   
1432    USE pegrid
1433   
1434    USE statistics
1435
1436    IMPLICIT NONE
1437
1438    INTEGER(iwp) ::  i                   !:
1439    INTEGER(iwp) ::  j                   !:
1440    INTEGER(iwp) ::  k                   !:
1441    INTEGER(iwp) ::  nt                  !:
1442    INTEGER(iwp) ::  omp_get_thread_num  !:
1443    INTEGER(iwp) ::  sr                  !:
1444    INTEGER(iwp) ::  tn                  !:
1445   
1446    LOGICAL ::  first  !:
1447   
1448    REAL(wp) ::  dptdz_threshold  !:
1449    REAL(wp) ::  fac              !:
1450    REAL(wp) ::  height           !:
1451    REAL(wp) ::  pts              !:
1452    REAL(wp) ::  sums_l_eper      !:
1453    REAL(wp) ::  sums_l_etot      !:
1454    REAL(wp) ::  s1               !:
1455    REAL(wp) ::  s2               !:
1456    REAL(wp) ::  s3               !:
1457    REAL(wp) ::  s4               !:
1458    REAL(wp) ::  s5               !:
1459    REAL(wp) ::  s6               !:
1460    REAL(wp) ::  s7               !:
1461    REAL(wp) ::  ust              !:
1462    REAL(wp) ::  ust2             !:
1463    REAL(wp) ::  u2               !:
1464    REAL(wp) ::  vst              !:
1465    REAL(wp) ::  vst2             !:
1466    REAL(wp) ::  v2               !:
1467    REAL(wp) ::  w2               !:
1468    REAL(wp) ::  z_i(2)           !:
1469
1470    REAL(wp) ::  dptdz(nzb+1:nzt+1)    !:
1471    REAL(wp) ::  sums_ll(nzb:nzt+1,2)  !:
1472
1473    CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
1474
1475!
1476!-- To be on the safe side, check whether flow_statistics has already been
1477!-- called once after the current time step
1478    IF ( flow_statistics_called )  THEN
1479
1480       message_string = 'flow_statistics is called two times within one ' // &
1481                        'timestep'
1482       CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
1483
1484    ENDIF
1485
1486    !$acc data copyin( hom ) create( sums, sums_l )
1487
1488!
1489!-- Compute statistics for each (sub-)region
1490    DO  sr = 0, statistic_regions
1491
1492!
1493!--    Initialize (local) summation array
1494       sums_l = 0.0_wp
1495
1496!
1497!--    Store sums that have been computed in other subroutines in summation
1498!--    array
1499       sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
1500!--    WARNING: next line still has to be adjusted for OpenMP
1501       sums_l(:,21,0) = sums_wsts_bc_l(:,sr)  ! heat flux from advec_s_bc
1502       sums_l(nzb+9,pr_palm,0)  = sums_divold_l(sr)  ! old divergence from pres
1503       sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
1504
1505!
1506!--    Copy the turbulent quantities, evaluated in the advection routines to
1507!--    the local array sums_l() for further computations
1508       IF ( ws_scheme_mom .AND. sr == 0 )  THEN
1509
1510!
1511!--       According to the Neumann bc for the horizontal velocity components,
1512!--       the corresponding fluxes has to satisfiy the same bc.
1513          IF ( ocean )  THEN
1514             sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
1515             sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)
1516          ENDIF
1517
1518          DO  i = 0, threads_per_task-1
1519!
1520!--          Swap the turbulent quantities evaluated in advec_ws.
1521             sums_l(:,13,i) = sums_wsus_ws_l(:,i)       ! w*u*
1522             sums_l(:,15,i) = sums_wsvs_ws_l(:,i)       ! w*v*
1523             sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2
1524             sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2
1525             sums_l(:,32,i) = sums_ws2_ws_l(:,i)        ! w*2
1526             sums_l(:,34,i) = sums_l(:,34,i) + 0.5_wp *                        &
1527                              ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +      &
1528                                sums_ws2_ws_l(:,i) )    ! e*
1529             DO  k = nzb, nzt
1530                sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5_wp * ( &
1531                                                      sums_us2_ws_l(k,i) +     &
1532                                                      sums_vs2_ws_l(k,i) +     &
1533                                                      sums_ws2_ws_l(k,i)     )
1534             ENDDO
1535          ENDDO
1536
1537       ENDIF
1538
1539       IF ( ws_scheme_sca .AND. sr == 0 )  THEN
1540
1541          DO  i = 0, threads_per_task-1
1542             sums_l(:,17,i) = sums_wspts_ws_l(:,i)      ! w*pt* from advec_s_ws
1543             IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
1544             IF ( humidity .OR. passive_scalar ) sums_l(:,49,i) =              &
1545                                                   sums_wsqs_ws_l(:,i) !w*q*
1546          ENDDO
1547
1548       ENDIF
1549!
1550!--    Horizontally averaged profiles of horizontal velocities and temperature.
1551!--    They must have been computed before, because they are already required
1552!--    for other horizontal averages.
1553       tn = 0
1554
1555       !$OMP PARALLEL PRIVATE( i, j, k, tn )
1556#if defined( __intel_openmp_bug )
1557       tn = omp_get_thread_num()
1558#else
1559!$     tn = omp_get_thread_num()
1560#endif
1561
1562       !$acc update device( sums_l )
1563
1564       !$OMP DO
1565       !$acc parallel loop gang present( pt, rflags_invers, rmask, sums_l, u, v ) create( s1, s2, s3 )
1566       DO  k = nzb, nzt+1
1567          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
1568          DO  i = nxl, nxr
1569             DO  j =  nys, nyn
1570!
1571!--             k+1 is used in rflags since rflags is set 0 at surface points
1572                s1 = s1 + u(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1573                s2 = s2 + v(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1574                s3 = s3 + pt(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1575             ENDDO
1576          ENDDO
1577          sums_l(k,1,tn) = s1
1578          sums_l(k,2,tn) = s2
1579          sums_l(k,4,tn) = s3
1580       ENDDO
1581       !$acc end parallel loop
1582
1583!
1584!--    Horizontally averaged profile of salinity
1585       IF ( ocean )  THEN
1586          !$OMP DO
1587          !$acc parallel loop gang present( rflags_invers, rmask, sums_l, sa ) create( s1 )
1588          DO  k = nzb, nzt+1
1589             !$acc loop vector collapse( 2 ) reduction( +: s1 )
1590             DO  i = nxl, nxr
1591                DO  j =  nys, nyn
1592                   s1 = s1 + sa(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1593                ENDDO
1594             ENDDO
1595             sums_l(k,23,tn) = s1
1596          ENDDO
1597          !$acc end parallel loop
1598       ENDIF
1599
1600!
1601!--    Horizontally averaged profiles of virtual potential temperature,
1602!--    total water content, specific humidity and liquid water potential
1603!--    temperature
1604       IF ( humidity )  THEN
1605
1606          !$OMP DO
1607          !$acc parallel loop gang present( q, rflags_invers, rmask, sums_l, vpt ) create( s1, s2 )
1608          DO  k = nzb, nzt+1
1609             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
1610             DO  i = nxl, nxr
1611                DO  j =  nys, nyn
1612                   s1 = s1 + q(k,j,i)   * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1613                   s2 = s2 + vpt(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1614                ENDDO
1615             ENDDO
1616             sums_l(k,41,tn) = s1
1617             sums_l(k,44,tn) = s2
1618          ENDDO
1619          !$acc end parallel loop
1620
1621          IF ( cloud_physics )  THEN
1622             !$OMP DO
1623             !$acc parallel loop gang present( pt, q, ql, rflags_invers, rmask, sums_l ) create( s1, s2 )
1624             DO  k = nzb, nzt+1
1625                !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
1626                DO  i = nxl, nxr
1627                   DO  j =  nys, nyn
1628                      s1 = s1 + ( q(k,j,i) - ql(k,j,i) ) * &
1629                                rmask(j,i,sr) * rflags_invers(j,i,k+1)
1630                      s2 = s2 + ( pt(k,j,i) + l_d_cp*pt_d_t(k) * ql(k,j,i) ) * &
1631                                rmask(j,i,sr) * rflags_invers(j,i,k+1)
1632                   ENDDO
1633                ENDDO
1634                sums_l(k,42,tn) = s1
1635                sums_l(k,43,tn) = s2
1636             ENDDO
1637             !$acc end parallel loop
1638          ENDIF
1639       ENDIF
1640
1641!
1642!--    Horizontally averaged profiles of passive scalar
1643       IF ( passive_scalar )  THEN
1644          !$OMP DO
1645          !$acc parallel loop gang present( q, rflags_invers, rmask, sums_l ) create( s1 )
1646          DO  k = nzb, nzt+1
1647             !$acc loop vector collapse( 2 ) reduction( +: s1 )
1648             DO  i = nxl, nxr
1649                DO  j =  nys, nyn
1650                   s1 = s1 + q(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1651                ENDDO
1652             ENDDO
1653             sums_l(k,41,tn) = s1
1654          ENDDO
1655          !$acc end parallel loop
1656       ENDIF
1657       !$OMP END PARALLEL
1658
1659!
1660!--    Summation of thread sums
1661       IF ( threads_per_task > 1 )  THEN
1662          DO  i = 1, threads_per_task-1
1663             !$acc parallel present( sums_l )
1664             sums_l(:,1,0) = sums_l(:,1,0) + sums_l(:,1,i)
1665             sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i)
1666             sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i)
1667             !$acc end parallel
1668             IF ( ocean )  THEN
1669                !$acc parallel present( sums_l )
1670                sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)
1671                !$acc end parallel
1672             ENDIF
1673             IF ( humidity )  THEN
1674                !$acc parallel present( sums_l )
1675                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
1676                sums_l(:,44,0) = sums_l(:,44,0) + sums_l(:,44,i)
1677                !$acc end parallel
1678                IF ( cloud_physics )  THEN
1679                   !$acc parallel present( sums_l )
1680                   sums_l(:,42,0) = sums_l(:,42,0) + sums_l(:,42,i)
1681                   sums_l(:,43,0) = sums_l(:,43,0) + sums_l(:,43,i)
1682                   !$acc end parallel
1683                ENDIF
1684             ENDIF
1685             IF ( passive_scalar )  THEN
1686                !$acc parallel present( sums_l )
1687                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
1688                !$acc end parallel
1689             ENDIF
1690          ENDDO
1691       ENDIF
1692
1693#if defined( __parallel )
1694!
1695!--    Compute total sum from local sums
1696       !$acc update host( sums_l )
1697       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1698       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL,  &
1699                           MPI_SUM, comm2d, ierr )
1700       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1701       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL,  &
1702                           MPI_SUM, comm2d, ierr )
1703       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1704       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL,  &
1705                           MPI_SUM, comm2d, ierr )
1706       IF ( ocean )  THEN
1707          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1708          CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb,       &
1709                              MPI_REAL, MPI_SUM, comm2d, ierr )
1710       ENDIF
1711       IF ( humidity ) THEN
1712          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1713          CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb,       &
1714                              MPI_REAL, MPI_SUM, comm2d, ierr )
1715          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1716          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb,       &
1717                              MPI_REAL, MPI_SUM, comm2d, ierr )
1718          IF ( cloud_physics ) THEN
1719             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1720             CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb,    &
1721                                 MPI_REAL, MPI_SUM, comm2d, ierr )
1722             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1723             CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb,    &
1724                                 MPI_REAL, MPI_SUM, comm2d, ierr )
1725          ENDIF
1726       ENDIF
1727
1728       IF ( passive_scalar )  THEN
1729          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1730          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb,       &
1731                              MPI_REAL, MPI_SUM, comm2d, ierr )
1732       ENDIF
1733       !$acc update device( sums )
1734#else
1735       !$acc parallel present( sums, sums_l )
1736       sums(:,1) = sums_l(:,1,0)
1737       sums(:,2) = sums_l(:,2,0)
1738       sums(:,4) = sums_l(:,4,0)
1739       !$acc end parallel
1740       IF ( ocean )  THEN
1741          !$acc parallel present( sums, sums_l )
1742          sums(:,23) = sums_l(:,23,0)
1743          !$acc end parallel
1744       ENDIF
1745       IF ( humidity )  THEN
1746          !$acc parallel present( sums, sums_l )
1747          sums(:,44) = sums_l(:,44,0)
1748          sums(:,41) = sums_l(:,41,0)
1749          !$acc end parallel
1750          IF ( cloud_physics )  THEN
1751             !$acc parallel present( sums, sums_l )
1752             sums(:,42) = sums_l(:,42,0)
1753             sums(:,43) = sums_l(:,43,0)
1754             !$acc end parallel
1755          ENDIF
1756       ENDIF
1757       IF ( passive_scalar )  THEN
1758          !$acc parallel present( sums, sums_l )
1759          sums(:,41) = sums_l(:,41,0)
1760          !$acc end parallel
1761       ENDIF
1762#endif
1763
1764!
1765!--    Final values are obtained by division by the total number of grid points
1766!--    used for summation. After that store profiles.
1767       !$acc parallel present( hom, ngp_2dh, ngp_2dh_s_inner, sums )
1768       sums(:,1) = sums(:,1) / ngp_2dh(sr)
1769       sums(:,2) = sums(:,2) / ngp_2dh(sr)
1770       sums(:,4) = sums(:,4) / ngp_2dh_s_inner(:,sr)
1771       hom(:,1,1,sr) = sums(:,1)             ! u
1772       hom(:,1,2,sr) = sums(:,2)             ! v
1773       hom(:,1,4,sr) = sums(:,4)             ! pt
1774       !$acc end parallel
1775
1776!
1777!--    Salinity
1778       IF ( ocean )  THEN
1779          !$acc parallel present( hom, ngp_2dh_s_inner, sums )
1780          sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr)
1781          hom(:,1,23,sr) = sums(:,23)             ! sa
1782          !$acc end parallel
1783       ENDIF
1784
1785!
1786!--    Humidity and cloud parameters
1787       IF ( humidity ) THEN
1788          !$acc parallel present( hom, ngp_2dh_s_inner, sums )
1789          sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr)
1790          sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
1791          hom(:,1,44,sr) = sums(:,44)                ! vpt
1792          hom(:,1,41,sr) = sums(:,41)                ! qv (q)
1793          !$acc end parallel
1794          IF ( cloud_physics ) THEN
1795             !$acc parallel present( hom, ngp_2dh_s_inner, sums )
1796             sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr)
1797             sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr)
1798             hom(:,1,42,sr) = sums(:,42)             ! qv
1799             hom(:,1,43,sr) = sums(:,43)             ! pt
1800             !$acc end parallel
1801          ENDIF
1802       ENDIF
1803
1804!
1805!--    Passive scalar
1806       IF ( passive_scalar )  THEN
1807          !$acc parallel present( hom, ngp_2dh_s_inner, sums )
1808          sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
1809          hom(:,1,41,sr) = sums(:,41)                ! s (q)
1810          !$acc end parallel
1811       ENDIF
1812
1813!
1814!--    Horizontally averaged profiles of the remaining prognostic variables,
1815!--    variances, the total and the perturbation energy (single values in last
1816!--    column of sums_l) and some diagnostic quantities.
1817!--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
1818!--    ----  speaking the following k-loop would have to be split up and
1819!--          rearranged according to the staggered grid.
1820!--          However, this implies no error since staggered velocity components
1821!--          are zero at the walls and inside buildings.
1822       tn = 0
1823#if defined( __intel_openmp_bug )
1824       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, &
1825       !$OMP                    tn, ust, ust2, u2, vst, vst2, v2, w2 )
1826       tn = omp_get_thread_num()
1827#else
1828       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2, w2 )
1829!$     tn = omp_get_thread_num()
1830#endif
1831       !$OMP DO
1832       !$acc parallel loop gang present( e, hom, kh, km, p, pt, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3, s4, s5, s6, s7 )
1833       DO  k = nzb, nzt+1
1834          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5, s6, s7 )
1835          DO  i = nxl, nxr
1836             DO  j =  nys, nyn
1837!
1838!--             Prognostic and diagnostic variables
1839                s1 = s1 + w(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1840                s2 = s2 + e(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1841                s3 = s3 + km(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1842                s4 = s4 + kh(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1843                s5 = s5 + p(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1844                s6 = s6 + ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr) * &
1845                          rflags_invers(j,i,k+1)
1846!
1847!--             Higher moments
1848!--             (Computation of the skewness of w further below)
1849                s7 = s7 + w(k,j,i)**3 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1850             ENDDO
1851          ENDDO
1852          sums_l(k,3,tn)  = s1
1853          sums_l(k,8,tn)  = s2
1854          sums_l(k,9,tn)  = s3
1855          sums_l(k,10,tn) = s4
1856          sums_l(k,40,tn) = s5
1857          sums_l(k,33,tn) = s6
1858          sums_l(k,38,tn) = s7
1859       ENDDO
1860       !$acc end parallel loop
1861
1862       IF ( humidity )  THEN
1863          !$OMP DO
1864          !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l ) create( s1 )
1865          DO  k = nzb, nzt+1
1866             !$acc loop vector collapse( 2 ) reduction( +: s1 )
1867             DO  i = nxl, nxr
1868                DO  j =  nys, nyn
1869                   s1 = s1 + ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr) * &
1870                             rflags_invers(j,i,k+1)
1871                ENDDO
1872             ENDDO
1873             sums_l(k,70,tn) = s1
1874          ENDDO
1875          !$acc end parallel loop
1876       ENDIF
1877
1878!
1879!--    Total and perturbation energy for the total domain (being
1880!--    collected in the last column of sums_l).
1881       !$OMP DO
1882       !$acc parallel loop collapse(3) present( rflags_invers, rmask, u, v, w ) reduction(+:s1)
1883       DO  i = nxl, nxr
1884          DO  j =  nys, nyn
1885             DO  k = nzb, nzt+1
1886                s1 = s1 + 0.5_wp *                                             & 
1887                          ( u(k,j,i)**2 + v(k,j,i)**2 + w(k,j,i)**2 ) *        &
1888                          rmask(j,i,sr) * rflags_invers(j,i,k+1)
1889             ENDDO
1890          ENDDO
1891       ENDDO
1892       !$acc end parallel loop
1893       !$acc parallel present( sums_l )
1894       sums_l(nzb+4,pr_palm,tn) = s1
1895       !$acc end parallel
1896
1897       !$OMP DO
1898       !$acc parallel present( rmask, sums_l, us, usws, vsws, ts ) create( s1, s2, s3, s4 )
1899       !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4 )
1900       DO  i = nxl, nxr
1901          DO  j =  nys, nyn
1902!
1903!--          2D-arrays (being collected in the last column of sums_l)
1904             s1 = s1 + us(j,i)   * rmask(j,i,sr)
1905             s2 = s2 + usws(j,i) * rmask(j,i,sr)
1906             s3 = s3 + vsws(j,i) * rmask(j,i,sr)
1907             s4 = s4 + ts(j,i)   * rmask(j,i,sr)
1908          ENDDO
1909       ENDDO
1910       sums_l(nzb,pr_palm,tn)   = s1
1911       sums_l(nzb+1,pr_palm,tn) = s2
1912       sums_l(nzb+2,pr_palm,tn) = s3
1913       sums_l(nzb+3,pr_palm,tn) = s4
1914       !$acc end parallel
1915
1916       IF ( humidity )  THEN
1917          !$acc parallel present( qs, rmask, sums_l ) create( s1 )
1918          !$acc loop vector collapse( 2 ) reduction( +: s1 )
1919          DO  i = nxl, nxr
1920             DO  j =  nys, nyn
1921                s1 = s1 + qs(j,i) * rmask(j,i,sr)
1922             ENDDO
1923          ENDDO
1924          sums_l(nzb+12,pr_palm,tn) = s1
1925          !$acc end parallel
1926       ENDIF
1927
1928!
1929!--    Computation of statistics when ws-scheme is not used. Else these
1930!--    quantities are evaluated in the advection routines.
1931       IF ( .NOT. ws_scheme_mom  .OR.  sr /= 0 )  THEN
1932
1933          !$OMP DO
1934          !$acc parallel loop gang present( u, v, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3, s4, ust2, vst2, w2 )
1935          DO  k = nzb, nzt+1
1936             !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4 )
1937             DO  i = nxl, nxr
1938                DO  j =  nys, nyn
1939                   ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
1940                   vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
1941                   w2   = w(k,j,i)**2
1942
1943                   s1 = s1 + ust2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1944                   s2 = s2 + vst2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1945                   s3 = s3 + w2   * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1946!
1947!--                Perturbation energy
1948                   s4 = s4 + 0.5_wp * ( ust2 + vst2 + w2 ) * rmask(j,i,sr) *   &
1949                             rflags_invers(j,i,k+1)
1950                ENDDO
1951             ENDDO
1952             sums_l(k,30,tn) = s1
1953             sums_l(k,31,tn) = s2
1954             sums_l(k,32,tn) = s3
1955             sums_l(k,34,tn) = s4
1956          ENDDO
1957          !$acc end parallel loop
1958!
1959!--       Total perturbation TKE
1960          !$OMP DO
1961          !$acc parallel present( sums_l ) create( s1 )
1962          !$acc loop reduction( +: s1 )
1963          DO  k = nzb, nzt+1
1964             s1 = s1 + sums_l(k,34,tn)
1965          ENDDO
1966          sums_l(nzb+5,pr_palm,tn) = s1
1967          !$acc end parallel
1968
1969       ENDIF
1970
1971!
1972!--    Horizontally averaged profiles of the vertical fluxes
1973
1974!
1975!--    Subgridscale fluxes.
1976!--    WARNING: If a Prandtl-layer is used (k=nzb for flat terrain), the fluxes
1977!--    -------  should be calculated there in a different way. This is done
1978!--             in the next loop further below, where results from this loop are
1979!--             overwritten. However, THIS WORKS IN CASE OF FLAT TERRAIN ONLY!
1980!--             The non-flat case still has to be handled.
1981!--    NOTE: for simplicity, nzb_s_inner is used below, although
1982!--    ----  strictly speaking the following k-loop would have to be
1983!--          split up according to the staggered grid.
1984!--          However, this implies no error since staggered velocity
1985!--          components are zero at the walls and inside buildings.
1986       !$OMP DO
1987       !$acc parallel loop gang present( ddzu, kh, km, pt, u, v, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3 )
1988       DO  k = nzb, nzt_diff
1989          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
1990          DO  i = nxl, nxr
1991             DO  j = nys, nyn
1992
1993!
1994!--             Momentum flux w"u"
1995                s1 = s1 - 0.25_wp * (                                          &
1996                               km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
1997                                                           ) * (               &
1998                                   ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
1999                                 + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
2000                                                               )               &
2001                               * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2002!
2003!--             Momentum flux w"v"
2004                s2 = s2 - 0.25_wp * (                                          &
2005                               km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
2006                                                           ) * (               &
2007                                   ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
2008                                 + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
2009                                                               )               &
2010                               * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2011!
2012!--             Heat flux w"pt"
2013                s3 = s3 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )                 &
2014                                 * ( pt(k+1,j,i) - pt(k,j,i) )                 &
2015                                 * ddzu(k+1) * rmask(j,i,sr)                   &
2016                                 * rflags_invers(j,i,k+1)
2017             ENDDO
2018          ENDDO
2019          sums_l(k,12,tn) = s1
2020          sums_l(k,14,tn) = s2
2021          sums_l(k,16,tn) = s3
2022       ENDDO
2023       !$acc end parallel loop
2024
2025!
2026!--    Salinity flux w"sa"
2027       IF ( ocean )  THEN
2028          !$acc parallel loop gang present( ddzu, kh, sa, rflags_invers, rmask, sums_l ) create( s1 )
2029          DO  k = nzb, nzt_diff
2030             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2031             DO  i = nxl, nxr
2032                DO  j = nys, nyn
2033                   s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
2034                                    * ( sa(k+1,j,i) - sa(k,j,i) )              &
2035                                    * ddzu(k+1) * rmask(j,i,sr)                & 
2036                                    * rflags_invers(j,i,k+1)
2037                ENDDO
2038             ENDDO
2039             sums_l(k,65,tn) = s1
2040          ENDDO
2041          !$acc end parallel loop
2042       ENDIF
2043
2044!
2045!--    Buoyancy flux, water flux (humidity flux) w"q"
2046       IF ( humidity ) THEN
2047
2048          !$acc parallel loop gang present( ddzu, kh, q, vpt, rflags_invers, rmask, sums_l ) create( s1, s2 )
2049          DO  k = nzb, nzt_diff
2050             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2051             DO  i = nxl, nxr
2052                DO  j = nys, nyn
2053                   s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
2054                                    * ( vpt(k+1,j,i) - vpt(k,j,i) )            &
2055                                    * ddzu(k+1) * rmask(j,i,sr)                &
2056                                    * rflags_invers(j,i,k+1)
2057                   s2 = s2 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
2058                                    * ( q(k+1,j,i) - q(k,j,i) )                &
2059                                    * ddzu(k+1) * rmask(j,i,sr)                &
2060                                    * rflags_invers(j,i,k+1)
2061                ENDDO
2062             ENDDO
2063             sums_l(k,45,tn) = s1
2064             sums_l(k,48,tn) = s2
2065          ENDDO
2066          !$acc end parallel loop
2067
2068          IF ( cloud_physics ) THEN
2069
2070             !$acc parallel loop gang present( ddzu, kh, q, ql, rflags_invers, rmask, sums_l ) create( s1 )
2071             DO  k = nzb, nzt_diff
2072                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2073                DO  i = nxl, nxr
2074                   DO  j = nys, nyn
2075                      s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )           &
2076                                       *  ( ( q(k+1,j,i) - ql(k+1,j,i) )       &
2077                                          - ( q(k,j,i) - ql(k,j,i) ) )         &
2078                                       * ddzu(k+1) * rmask(j,i,sr)             & 
2079                                       * rflags_invers(j,i,k+1)
2080                   ENDDO
2081                ENDDO
2082                sums_l(k,51,tn) = s1
2083             ENDDO
2084             !$acc end parallel loop
2085
2086          ENDIF
2087
2088       ENDIF
2089!
2090!--    Passive scalar flux
2091       IF ( passive_scalar )  THEN
2092
2093          !$acc parallel loop gang present( ddzu, kh, q, rflags_invers, rmask, sums_l ) create( s1 )
2094          DO  k = nzb, nzt_diff
2095             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2096             DO  i = nxl, nxr
2097                DO  j = nys, nyn
2098                   s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
2099                                    * ( q(k+1,j,i) - q(k,j,i) )                &
2100                                    * ddzu(k+1) * rmask(j,i,sr)                &
2101                                    * rflags_invers(j,i,k+1)
2102                ENDDO
2103             ENDDO
2104             sums_l(k,48,tn) = s1
2105          ENDDO
2106          !$acc end parallel loop
2107
2108       ENDIF
2109
2110       IF ( use_surface_fluxes )  THEN
2111
2112          !$OMP DO
2113          !$acc parallel present( rmask, shf, sums_l, usws, vsws ) create( s1, s2, s3, s4, s5 )
2114          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5 )
2115          DO  i = nxl, nxr
2116             DO  j =  nys, nyn
2117!
2118!--             Subgridscale fluxes in the Prandtl layer
2119                s1 = s1 + usws(j,i) * rmask(j,i,sr)     ! w"u"
2120                s2 = s2 + vsws(j,i) * rmask(j,i,sr)     ! w"v"
2121                s3 = s3 + shf(j,i)  * rmask(j,i,sr)     ! w"pt"
2122                s4 = s4 + 0.0_wp * rmask(j,i,sr)        ! u"pt"
2123                s5 = s5 + 0.0_wp * rmask(j,i,sr)        ! v"pt"
2124             ENDDO
2125          ENDDO
2126          sums_l(nzb,12,tn) = s1
2127          sums_l(nzb,14,tn) = s2
2128          sums_l(nzb,16,tn) = s3
2129          sums_l(nzb,58,tn) = s4
2130          sums_l(nzb,61,tn) = s5
2131          !$acc end parallel
2132
2133          IF ( ocean )  THEN
2134
2135             !$OMP DO
2136             !$acc parallel present( rmask, saswsb, sums_l ) create( s1 )
2137             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2138             DO  i = nxl, nxr
2139                DO  j =  nys, nyn
2140                   s1 = s1 + saswsb(j,i) * rmask(j,i,sr)  ! w"sa"
2141                ENDDO
2142             ENDDO
2143             sums_l(nzb,65,tn) = s1
2144             !$acc end parallel
2145
2146          ENDIF
2147
2148          IF ( humidity )  THEN
2149
2150             !$OMP DO
2151             !$acc parallel present( pt, q, qsws, rmask, shf, sums_l ) create( s1, s2 )
2152             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2153             DO  i = nxl, nxr
2154                DO  j =  nys, nyn
2155                   s1 = s1 + qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
2156                   s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) * shf(j,i)    &
2157                               + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )
2158                ENDDO
2159             ENDDO
2160             sums_l(nzb,48,tn) = s1
2161             sums_l(nzb,45,tn) = s2
2162             !$acc end parallel
2163
2164             IF ( cloud_droplets )  THEN
2165
2166                !$OMP DO
2167                !$acc parallel present( pt, q, ql, qsws, rmask, shf, sums_l ) create( s1 )
2168                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2169                DO  i = nxl, nxr
2170                   DO  j =  nys, nyn
2171                      s1 = s1 + ( ( 1.0_wp +                                   &
2172                                    0.61_wp * q(nzb,j,i) - ql(nzb,j,i) ) *     &
2173                                  shf(j,i) + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )
2174                   ENDDO
2175                ENDDO
2176                sums_l(nzb,45,tn) = s1
2177                !$acc end parallel
2178
2179             ENDIF
2180
2181             IF ( cloud_physics )  THEN
2182
2183                !$OMP DO
2184                !$acc parallel present( qsws, rmask, sums_l ) create( s1 )
2185                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2186                DO  i = nxl, nxr
2187                   DO  j =  nys, nyn
2188!
2189!--                   Formula does not work if ql(nzb) /= 0.0
2190                      s1 = s1 + qsws(j,i) * rmask(j,i,sr)   ! w"q" (w"qv")
2191                   ENDDO
2192                ENDDO
2193                sums_l(nzb,51,tn) = s1
2194                !$acc end parallel
2195
2196             ENDIF
2197
2198          ENDIF
2199
2200          IF ( passive_scalar )  THEN
2201
2202             !$OMP DO
2203             !$acc parallel present( qsws, rmask, sums_l ) create( s1 )
2204             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2205             DO  i = nxl, nxr
2206                DO  j =  nys, nyn
2207                   s1 = s1 + qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
2208                ENDDO
2209             ENDDO
2210             sums_l(nzb,48,tn) = s1
2211             !$acc end parallel
2212
2213          ENDIF
2214
2215       ENDIF
2216
2217!
2218!--    Subgridscale fluxes at the top surface
2219       IF ( use_top_fluxes )  THEN
2220
2221          !$OMP DO
2222          !$acc parallel present( rmask, sums_l, tswst, uswst, vswst ) create( s1, s2, s3, s4, s5 )
2223          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5 )
2224          DO  i = nxl, nxr
2225             DO  j =  nys, nyn
2226                s1 = s1 + uswst(j,i) * rmask(j,i,sr)    ! w"u"
2227                s2 = s2 + vswst(j,i) * rmask(j,i,sr)    ! w"v"
2228                s3 = s3 + tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
2229                s4 = s4 + 0.0_wp * rmask(j,i,sr)        ! u"pt"
2230                s5 = s5 + 0.0_wp * rmask(j,i,sr)        ! v"pt"
2231             ENDDO
2232          ENDDO
2233          sums_l(nzt:nzt+1,12,tn) = s1
2234          sums_l(nzt:nzt+1,14,tn) = s2
2235          sums_l(nzt:nzt+1,16,tn) = s3
2236          sums_l(nzt:nzt+1,58,tn) = s4
2237          sums_l(nzt:nzt+1,61,tn) = s5
2238          !$acc end parallel
2239
2240          IF ( ocean )  THEN
2241
2242             !$OMP DO
2243             !$acc parallel present( rmask, saswst, sums_l ) create( s1 )
2244             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2245             DO  i = nxl, nxr
2246                DO  j =  nys, nyn
2247                   s1 = s1 + saswst(j,i) * rmask(j,i,sr)  ! w"sa"
2248                ENDDO
2249             ENDDO
2250             sums_l(nzt,65,tn) = s1
2251             !$acc end parallel
2252
2253          ENDIF
2254
2255          IF ( humidity )  THEN
2256
2257             !$OMP DO
2258             !$acc parallel present( pt, q, qswst, rmask, tswst, sums_l ) create( s1, s2 )
2259             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2260             DO  i = nxl, nxr
2261                DO  j =  nys, nyn
2262                   s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
2263                   s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) * tswst(j,i) +&
2264                                 0.61_wp * pt(nzt,j,i) * qswst(j,i) )
2265                ENDDO
2266             ENDDO
2267             sums_l(nzt,48,tn) = s1
2268             sums_l(nzt,45,tn) = s2
2269             !$acc end parallel
2270
2271             IF ( cloud_droplets )  THEN
2272
2273                !$OMP DO
2274                !$acc parallel present( pt, q, ql, qswst, rmask, tswst, sums_l ) create( s1 )
2275                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2276                DO  i = nxl, nxr
2277                   DO  j =  nys, nyn
2278                      s1 = s1 + ( ( 1.0_wp +                                   &
2279                                    0.61_wp * q(nzt,j,i) - ql(nzt,j,i) ) *     &
2280                                  tswst(j,i) +                                 &
2281                                  0.61_wp * pt(nzt,j,i) * qswst(j,i) )
2282                   ENDDO
2283                ENDDO
2284                sums_l(nzt,45,tn) = s1
2285                !$acc end parallel
2286
2287             ENDIF
2288
2289             IF ( cloud_physics )  THEN
2290
2291                !$OMP DO
2292                !$acc parallel present( qswst, rmask, sums_l ) create( s1 )
2293                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2294                DO  i = nxl, nxr
2295                   DO  j =  nys, nyn
2296!
2297!--                   Formula does not work if ql(nzb) /= 0.0
2298                      s1 = s1 + qswst(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
2299                   ENDDO
2300                ENDDO
2301                sums_l(nzt,51,tn) = s1
2302                !$acc end parallel
2303
2304             ENDIF
2305
2306          ENDIF
2307
2308          IF ( passive_scalar )  THEN
2309
2310             !$OMP DO
2311             !$acc parallel present( qswst, rmask, sums_l ) create( s1 )
2312             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2313             DO  i = nxl, nxr
2314                DO  j =  nys, nyn
2315                   s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
2316                ENDDO
2317             ENDDO
2318             sums_l(nzt,48,tn) = s1
2319             !$acc end parallel
2320
2321          ENDIF
2322
2323       ENDIF
2324
2325!
2326!--    Resolved fluxes (can be computed for all horizontal points)
2327!--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
2328!--    ----  speaking the following k-loop would have to be split up and
2329!--          rearranged according to the staggered grid.
2330       !$acc parallel loop gang present( hom, pt, rflags_invers, rmask, sums_l, u, v, w ) create( s1, s2, s3 )
2331       DO  k = nzb, nzt_diff
2332          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
2333          DO  i = nxl, nxr
2334             DO  j = nys, nyn
2335                ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) + &
2336                                 u(k+1,j,i) - hom(k+1,1,1,sr) )
2337                vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) + &
2338                                 v(k+1,j,i) - hom(k+1,1,2,sr) )
2339                pts = 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) + &
2340                                 pt(k+1,j,i) - hom(k+1,1,4,sr) )
2341!
2342!--             Higher moments
2343                s1 = s1 + pts * w(k,j,i)**2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2344                s2 = s2 + pts**2 * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2345!
2346!--             Energy flux w*e* (has to be adjusted?)
2347                s3 = s3 + w(k,j,i) * 0.5_wp * ( ust**2 + vst**2 + w(k,j,i)**2 )&
2348                                   * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2349             ENDDO
2350          ENDDO
2351          sums_l(k,35,tn) = s1
2352          sums_l(k,36,tn) = s2
2353          sums_l(k,37,tn) = s3
2354       ENDDO
2355       !$acc end parallel loop
2356
2357!
2358!--    Salinity flux and density (density does not belong to here,
2359!--    but so far there is no other suitable place to calculate)
2360       IF ( ocean )  THEN
2361
2362          IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
2363
2364             !$acc parallel loop gang present( hom, rflags_invers, rmask, sa, sums_l, w ) create( s1 )
2365             DO  k = nzb, nzt_diff
2366                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2367                DO  i = nxl, nxr
2368                   DO  j = nys, nyn
2369                      s1 = s1 + 0.5_wp * ( sa(k,j,i)   - hom(k,1,23,sr) +      &
2370                                           sa(k+1,j,i) - hom(k+1,1,23,sr) )    &
2371                                       * w(k,j,i) * rmask(j,i,sr)              & 
2372                                       * rflags_invers(j,i,k+1)
2373                   ENDDO
2374                ENDDO
2375                sums_l(k,66,tn) = s1
2376             ENDDO
2377             !$acc end parallel loop
2378
2379          ENDIF
2380
2381          !$acc parallel loop gang present( rflags_invers, rho, prho, rmask, sums_l ) create( s1, s2 )
2382          DO  k = nzb, nzt_diff
2383             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2384             DO  i = nxl, nxr
2385                DO  j = nys, nyn
2386                   s1 = s1 + rho(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2387                   s2 = s2 + prho(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2388                ENDDO
2389             ENDDO
2390             sums_l(k,64,tn) = s1
2391             sums_l(k,71,tn) = s2
2392          ENDDO
2393          !$acc end parallel loop
2394
2395       ENDIF
2396
2397!
2398!--    Buoyancy flux, water flux, humidity flux, liquid water
2399!--    content, rain drop concentration and rain water content
2400       IF ( humidity )  THEN
2401
2402          IF ( cloud_physics  .OR.  cloud_droplets )  THEN
2403
2404             !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, vpt, w ) create( s1 )
2405             DO  k = nzb, nzt_diff
2406                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2407                DO  i = nxl, nxr
2408                   DO  j = nys, nyn
2409                      s1 = s1 + 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +     &
2410                                           vpt(k+1,j,i) - hom(k+1,1,44,sr) ) * &
2411                                         w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2412                   ENDDO
2413                ENDDO
2414                sums_l(k,46,tn) = s1
2415             ENDDO
2416             !$acc end parallel loop
2417
2418             IF ( .NOT. cloud_droplets )  THEN
2419
2420                !$acc parallel loop gang present( hom, q, ql, rflags_invers, rmask, sums_l, w ) create( s1 )
2421                DO  k = nzb, nzt_diff
2422                   !$acc loop vector collapse( 2 ) reduction( +: s1 )
2423                   DO  i = nxl, nxr
2424                      DO  j = nys, nyn
2425                         s1 = s1 + 0.5_wp * ( ( q(k,j,i)   - ql(k,j,i)   ) - hom(k,1,42,sr) +   &
2426                                              ( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) ) &
2427                                          * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2428                      ENDDO
2429                   ENDDO
2430                   sums_l(k,52,tn) = s1
2431                ENDDO
2432                !$acc end parallel loop
2433
2434                IF ( icloud_scheme == 0  )  THEN
2435
2436                   !$acc parallel loop gang present( qc, ql, rflags_invers, rmask, sums_l ) create( s1, s2 )
2437                   DO  k = nzb, nzt_diff
2438                      !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2439                      DO  i = nxl, nxr
2440                         DO  j = nys, nyn
2441                            s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2442                            s2 = s2 + qc(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2443                         ENDDO
2444                      ENDDO
2445                      sums_l(k,54,tn) = s1
2446                      sums_l(k,75,tn) = s2
2447                   ENDDO
2448                   !$acc end parallel loop
2449
2450                   IF ( precipitation )  THEN
2451
2452                      !$acc parallel loop gang present( nr, qr, prr, rflags_invers, rmask, sums_l ) create( s1, s2, s3 )
2453                      DO  k = nzb, nzt_diff
2454                         !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
2455                         DO  i = nxl, nxr
2456                            DO  j = nys, nyn
2457                               s1 = s1 + nr(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2458                               s2 = s2 + qr(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2459                               s3 = s3 + prr(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2460                            ENDDO
2461                         ENDDO
2462                         sums_l(k,73,tn) = s1
2463                         sums_l(k,74,tn) = s2
2464                         sums_l(k,76,tn) = s3
2465                      ENDDO
2466                      !$acc end parallel loop
2467
2468                   ENDIF
2469
2470                ELSE
2471
2472                   !$acc parallel loop gang present( ql, rflags_invers, rmask, sums_l ) create( s1 )
2473                   DO  k = nzb, nzt_diff
2474                      !$acc loop vector collapse( 2 ) reduction( +: s1 )
2475                      DO  i = nxl, nxr
2476                         DO  j = nys, nyn
2477                            s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2478                         ENDDO
2479                      ENDDO
2480                      sums_l(k,54,tn) = s1
2481                   ENDDO
2482                   !$acc end parallel loop
2483
2484                ENDIF
2485
2486             ELSE
2487
2488                !$acc parallel loop gang present( ql, rflags_invers, rmask, sums_l ) create( s1 )
2489                DO  k = nzb, nzt_diff
2490                   !$acc loop vector collapse( 2 ) reduction( +: s1 )
2491                   DO  i = nxl, nxr
2492                      DO  j = nys, nyn
2493                         s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2494                      ENDDO
2495                   ENDDO
2496                   sums_l(k,54,tn) = s1
2497                ENDDO
2498                !$acc end parallel loop
2499
2500             ENDIF
2501
2502          ELSE
2503
2504             IF( .NOT. ws_scheme_sca  .OR.  sr /= 0 )  THEN
2505
2506                !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, vpt, w ) create( s1 )
2507                DO  k = nzb, nzt_diff
2508                   !$acc loop vector collapse( 2 ) reduction( +: s1 )
2509                   DO  i = nxl, nxr
2510                      DO  j = nys, nyn
2511                         s1 = s1 + 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +   &
2512                                              vpt(k+1,j,i) - hom(k+1,1,44,sr) ) &
2513                                          * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2514                      ENDDO
2515                   ENDDO
2516                   sums_l(k,46,tn) = s1
2517                ENDDO
2518                !$acc end parallel loop
2519
2520             ELSEIF ( ws_scheme_sca  .AND.  sr == 0 )  THEN
2521
2522                !$acc parallel loop present( hom, sums_l )
2523                DO  k = nzb, nzt_diff
2524                   sums_l(k,46,tn) = ( 1.0_wp + 0.61_wp * hom(k,1,41,sr) ) * sums_l(k,17,tn) + &
2525                                                0.61_wp * hom(k,1,4,sr) * sums_l(k,49,tn)
2526                ENDDO
2527                !$acc end parallel loop
2528
2529             ENDIF
2530
2531          ENDIF
2532
2533       ENDIF
2534!
2535!--    Passive scalar flux
2536       IF ( passive_scalar  .AND.  ( .NOT. ws_scheme_sca  .OR.  sr /= 0 ) )  THEN
2537
2538          !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l, w ) create( s1 )
2539          DO  k = nzb, nzt_diff
2540             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2541             DO  i = nxl, nxr
2542                DO  j = nys, nyn
2543                   s1 = s1 + 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +          &
2544                                        q(k+1,j,i) - hom(k+1,1,41,sr) )        &
2545                                    * w(k,j,i) * rmask(j,i,sr)                 &
2546                                    * rflags_invers(j,i,k+1)
2547                ENDDO
2548             ENDDO
2549             sums_l(k,49,tn) = s1
2550          ENDDO
2551          !$acc end parallel loop
2552
2553       ENDIF
2554
2555!
2556!--    For speed optimization fluxes which have been computed in part directly
2557!--    inside the WS advection routines are treated seperatly
2558!--    Momentum fluxes first:
2559       IF ( .NOT. ws_scheme_mom  .OR.  sr /= 0  )  THEN
2560
2561          !$OMP DO
2562          !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, u, v, w ) create( s1, s2 )
2563          DO  k = nzb, nzt_diff
2564             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2565             DO  i = nxl, nxr
2566                DO  j = nys, nyn
2567                   ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +               &
2568                                    u(k+1,j,i) - hom(k+1,1,1,sr) )
2569                   vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +               &
2570                                    v(k+1,j,i) - hom(k+1,1,2,sr) )
2571!
2572!--                Momentum flux w*u*
2573                   s1 = s1 + 0.5_wp * ( w(k,j,i-1) + w(k,j,i) )                &
2574                                    * ust * rmask(j,i,sr)                      &
2575                                    * rflags_invers(j,i,k+1)
2576!
2577!--                Momentum flux w*v*
2578                   s2 = s2 + 0.5_wp * ( w(k,j-1,i) + w(k,j,i) )                &
2579                                    * vst * rmask(j,i,sr)                      & 
2580                                    * rflags_invers(j,i,k+1)
2581                ENDDO
2582             ENDDO
2583             sums_l(k,13,tn) = s1
2584             sums_l(k,15,tn) = s1
2585          ENDDO
2586          !$acc end parallel loop
2587
2588       ENDIF
2589
2590       IF ( .NOT. ws_scheme_sca  .OR.  sr /= 0 )  THEN
2591
2592          !$OMP DO
2593          !$acc parallel loop gang present( hom, pt, rflags_invers, rmask, sums_l, w ) create( s1 )
2594          DO  k = nzb, nzt_diff
2595             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2596             DO  i = nxl, nxr
2597                DO  j = nys, nyn
2598!
2599!--                Vertical heat flux
2600                   s1 = s1 + 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) +          &
2601                                        pt(k+1,j,i) - hom(k+1,1,4,sr) )        &
2602                                    * w(k,j,i) * rmask(j,i,sr)                 & 
2603                                    * rflags_invers(j,i,k+1)
2604                ENDDO
2605             ENDDO
2606             sums_l(k,17,tn) = s1
2607          ENDDO
2608          !$acc end parallel loop
2609
2610          IF ( humidity )  THEN
2611
2612             !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l, w ) create( s1 )
2613             DO  k = nzb, nzt_diff
2614                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2615                DO  i = nxl, nxr
2616                   DO  j = nys, nyn
2617                      s1 = s1 + 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +       &
2618                                           q(k+1,j,i) - hom(k+1,1,41,sr) )     &
2619                                       * w(k,j,i) * rmask(j,i,sr)              &
2620                                       * rflags_invers(j,i,k+1)
2621                   ENDDO
2622                ENDDO
2623                sums_l(k,49,tn) = s1
2624             ENDDO
2625             !$acc end parallel loop
2626
2627          ENDIF
2628
2629       ENDIF
2630
2631
2632!
2633!--    Density at top follows Neumann condition
2634       IF ( ocean )  THEN
2635          !$acc parallel present( sums_l )
2636          sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
2637          sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
2638          !$acc end parallel
2639       ENDIF
2640
2641!
2642!--    Divergence of vertical flux of resolved scale energy and pressure
2643!--    fluctuations as well as flux of pressure fluctuation itself (68).
2644!--    First calculate the products, then the divergence.
2645!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
2646       IF ( hom(nzb+1,2,55,0) /= 0.0_wp  .OR.  hom(nzb+1,2,68,0) /= 0.0_wp )  THEN
2647
2648          STOP '+++ openACC porting for vertical flux div of resolved scale TKE in flow_statistics is still missing'
2649          sums_ll = 0.0_wp  ! local array
2650
2651          !$OMP DO
2652          DO  i = nxl, nxr
2653             DO  j = nys, nyn
2654                DO  k = nzb_s_inner(j,i)+1, nzt
2655
2656                   sums_ll(k,1) = sums_ll(k,1) + 0.5_wp * w(k,j,i) * (         &
2657                  ( 0.25_wp * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1)    &
2658                              - 0.5_wp * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) )   &
2659                              ) )**2                                           &
2660                + ( 0.25_wp * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i)    &
2661                              - 0.5_wp * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) )   &
2662                              ) )**2                                           &
2663                + w(k,j,i)**2                                        )
2664
2665                   sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i)             &
2666                                               * ( p(k,j,i) + p(k+1,j,i) )
2667
2668                ENDDO
2669             ENDDO
2670          ENDDO
2671          sums_ll(0,1)     = 0.0_wp    ! because w is zero at the bottom
2672          sums_ll(nzt+1,1) = 0.0_wp
2673          sums_ll(0,2)     = 0.0_wp
2674          sums_ll(nzt+1,2) = 0.0_wp
2675
2676          DO  k = nzb+1, nzt
2677             sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
2678             sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
2679             sums_l(k,68,tn) = sums_ll(k,2)
2680          ENDDO
2681          sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
2682          sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
2683          sums_l(nzb,68,tn) = 0.0_wp    ! because w* = 0 at nzb
2684
2685       ENDIF
2686
2687!
2688!--    Divergence of vertical flux of SGS TKE and the flux itself (69)
2689       IF ( hom(nzb+1,2,57,0) /= 0.0_wp  .OR.  hom(nzb+1,2,69,0) /= 0.0_wp )  THEN
2690
2691          STOP '+++ openACC porting for vertical flux div of SGS TKE in flow_statistics is still missing'
2692          !$OMP DO
2693          DO  i = nxl, nxr
2694             DO  j = nys, nyn
2695                DO  k = nzb_s_inner(j,i)+1, nzt
2696
2697                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5_wp * (              &
2698                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
2699                 - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
2700                                                                ) * ddzw(k)
2701
2702                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5_wp * (              &
2703                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
2704                                                                )
2705
2706                ENDDO
2707             ENDDO
2708          ENDDO
2709          sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
2710          sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
2711
2712       ENDIF
2713
2714!
2715!--    Horizontal heat fluxes (subgrid, resolved, total).
2716!--    Do it only, if profiles shall be plotted.
2717       IF ( hom(nzb+1,2,58,0) /= 0.0_wp ) THEN
2718
2719          STOP '+++ openACC porting for horizontal flux calculation in flow_statistics is still missing'
2720          !$OMP DO
2721          DO  i = nxl, nxr
2722             DO  j = nys, nyn
2723                DO  k = nzb_s_inner(j,i)+1, nzt
2724!
2725!--                Subgrid horizontal heat fluxes u"pt", v"pt"
2726                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5_wp *                &
2727                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
2728                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
2729                                                 * ddx * rmask(j,i,sr)
2730                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
2731                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
2732                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
2733                                                 * ddy * rmask(j,i,sr)
2734!
2735!--                Resolved horizontal heat fluxes u*pt*, v*pt*
2736                   sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
2737                                     ( u(k,j,i) - hom(k,1,1,sr) ) * 0.5_wp *   &
2738                                     ( pt(k,j,i-1) - hom(k,1,4,sr) +           &
2739                                       pt(k,j,i)   - hom(k,1,4,sr) )
2740                   pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
2741                                    pt(k,j,i)   - hom(k,1,4,sr) )
2742                   sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
2743                                     ( v(k,j,i) - hom(k,1,2,sr) ) * 0.5_wp *   &
2744                                     ( pt(k,j-1,i) - hom(k,1,4,sr) +           &
2745                                       pt(k,j,i)   - hom(k,1,4,sr) )
2746                ENDDO
2747             ENDDO
2748          ENDDO
2749!
2750!--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
2751          sums_l(nzb,58,tn) = 0.0_wp
2752          sums_l(nzb,59,tn) = 0.0_wp
2753          sums_l(nzb,60,tn) = 0.0_wp
2754          sums_l(nzb,61,tn) = 0.0_wp
2755          sums_l(nzb,62,tn) = 0.0_wp
2756          sums_l(nzb,63,tn) = 0.0_wp
2757
2758       ENDIF
2759
2760!
2761!--    Collect current large scale advection and subsidence tendencies for
2762!--    data output
2763       IF ( large_scale_forcing )  THEN
2764!
2765!--       Interpolation in time of LSF_DATA
2766          nt = 1
2767          DO WHILE ( simulated_time - dt_3d > time_vert(nt) )
2768             nt = nt + 1
2769          ENDDO
2770          IF ( simulated_time - dt_3d /= time_vert(nt) )  THEN
2771            nt = nt - 1
2772          ENDIF
2773
2774          fac = ( simulated_time - dt_3d - time_vert(nt) )                     &
2775                / ( time_vert(nt+1)-time_vert(nt) )
2776
2777
2778          DO  k = nzb, nzt
2779             sums_ls_l(k,0) = td_lsa_lpt(k,nt)                                 &
2780                              + fac * ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )
2781             sums_ls_l(k,1) = td_lsa_q(k,nt)                                   &
2782                              + fac * ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) )
2783          ENDDO
2784
2785          sums_ls_l(nzt+1,0) = sums_ls_l(nzt,0)
2786          sums_ls_l(nzt+1,1) = sums_ls_l(nzt,1)
2787
2788          IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
2789
2790             DO  k = nzb, nzt
2791                sums_ls_l(k,2) = td_sub_lpt(k,nt) + fac *                      &
2792                                 ( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )
2793                sums_ls_l(k,3) = td_sub_q(k,nt) + fac *                        &
2794                                 ( td_sub_q(k,nt+1) - td_sub_q(k,nt) )
2795             ENDDO
2796
2797             sums_ls_l(nzt+1,2) = sums_ls_l(nzt,2)
2798             sums_ls_l(nzt+1,3) = sums_ls_l(nzt,3)
2799
2800          ENDIF
2801
2802       ENDIF
2803
2804!
2805!--    Calculate the user-defined profiles
2806       CALL user_statistics( 'profiles', sr, tn )
2807       !$OMP END PARALLEL
2808
2809!
2810!--    Summation of thread sums
2811       IF ( threads_per_task > 1 )  THEN
2812          STOP '+++ openACC porting for threads_per_task > 1 in flow_statistics is still missing'
2813          DO  i = 1, threads_per_task-1
2814             sums_l(:,3,0)          = sums_l(:,3,0) + sums_l(:,3,i)
2815             sums_l(:,4:40,0)       = sums_l(:,4:40,0) + sums_l(:,4:40,i)
2816             sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
2817                                      sums_l(:,45:pr_palm,i)
2818             IF ( max_pr_user > 0 )  THEN
2819                sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
2820                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
2821                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
2822             ENDIF
2823          ENDDO
2824       ENDIF
2825
2826       !$acc update host( hom, sums, sums_l )
2827
2828#if defined( __parallel )
2829
2830!
2831!--    Compute total sum from local sums
2832       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2833       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, &
2834                           MPI_SUM, comm2d, ierr )
2835       IF ( large_scale_forcing )  THEN
2836          CALL MPI_ALLREDUCE( sums_ls_l(nzb,2), sums(nzb,83), ngp_sums_ls,     &
2837                              MPI_REAL, MPI_SUM, comm2d, ierr )
2838       ENDIF
2839#else
2840       sums = sums_l(:,:,0)
2841       IF ( large_scale_forcing )  THEN
2842          sums(:,81:88) = sums_ls_l
2843       ENDIF
2844#endif
2845
2846!
2847!--    Final values are obtained by division by the total number of grid points
2848!--    used for summation. After that store profiles.
2849!--    Profiles:
2850       DO  k = nzb, nzt+1
2851          sums(k,3)               = sums(k,3)           / ngp_2dh(sr)
2852          sums(k,8:11)            = sums(k,8:11)        / ngp_2dh_s_inner(k,sr)
2853          sums(k,12:22)           = sums(k,12:22)       / ngp_2dh(sr)
2854          sums(k,23:29)           = sums(k,23:29)       / ngp_2dh_s_inner(k,sr)
2855          sums(k,30:32)           = sums(k,30:32)       / ngp_2dh(sr)
2856          sums(k,33:34)           = sums(k,33:34)       / ngp_2dh_s_inner(k,sr)
2857          sums(k,35:39)           = sums(k,35:39)       / ngp_2dh(sr)
2858          sums(k,40)              = sums(k,40)          / ngp_2dh_s_inner(k,sr)
2859          sums(k,45:53)           = sums(k,45:53)       / ngp_2dh(sr)
2860          sums(k,54)              = sums(k,54)          / ngp_2dh_s_inner(k,sr)
2861          sums(k,55:63)           = sums(k,55:63)       / ngp_2dh(sr)
2862          sums(k,64)              = sums(k,64)          / ngp_2dh_s_inner(k,sr)
2863          sums(k,70:80)           = sums(k,70:80)       / ngp_2dh_s_inner(k,sr)
2864          sums(k,81:88)           = sums(k,81:88)       / ngp_2dh(sr)
2865          sums(k,89:pr_palm-2)    = sums(k,89:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
2866       ENDDO
2867
2868!--    Upstream-parts
2869       sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
2870!--    u* and so on
2871!--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
2872!--    size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
2873!--    above the topography, they are being divided by ngp_2dh(sr)
2874       sums(nzb:nzb+3,pr_palm)    = sums(nzb:nzb+3,pr_palm)    / &
2875                                    ngp_2dh(sr)
2876       sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)       / &    ! qs
2877                                    ngp_2dh(sr)
2878!--    eges, e*
2879       sums(nzb+4:nzb+5,pr_palm)  = sums(nzb+4:nzb+5,pr_palm)  / &
2880                                    ngp_3d(sr)
2881!--    Old and new divergence
2882       sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
2883                                    ngp_3d_inner(sr)
2884
2885!--    User-defined profiles
2886       IF ( max_pr_user > 0 )  THEN
2887          DO  k = nzb, nzt+1
2888             sums(k,pr_palm+1:pr_palm+max_pr_user) = &
2889                                    sums(k,pr_palm+1:pr_palm+max_pr_user) / &
2890                                    ngp_2dh_s_inner(k,sr)
2891          ENDDO
2892       ENDIF
2893
2894!
2895!--    Collect horizontal average in hom.
2896!--    Compute deduced averages (e.g. total heat flux)
2897       hom(:,1,3,sr)  = sums(:,3)      ! w
2898       hom(:,1,8,sr)  = sums(:,8)      ! e     profiles 5-7 are initial profiles
2899       hom(:,1,9,sr)  = sums(:,9)      ! km
2900       hom(:,1,10,sr) = sums(:,10)     ! kh
2901       hom(:,1,11,sr) = sums(:,11)     ! l
2902       hom(:,1,12,sr) = sums(:,12)     ! w"u"
2903       hom(:,1,13,sr) = sums(:,13)     ! w*u*
2904       hom(:,1,14,sr) = sums(:,14)     ! w"v"
2905       hom(:,1,15,sr) = sums(:,15)     ! w*v*
2906       hom(:,1,16,sr) = sums(:,16)     ! w"pt"
2907       hom(:,1,17,sr) = sums(:,17)     ! w*pt*
2908       hom(:,1,18,sr) = sums(:,16) + sums(:,17)    ! wpt
2909       hom(:,1,19,sr) = sums(:,12) + sums(:,13)    ! wu
2910       hom(:,1,20,sr) = sums(:,14) + sums(:,15)    ! wv
2911       hom(:,1,21,sr) = sums(:,21)     ! w*pt*BC
2912       hom(:,1,22,sr) = sums(:,16) + sums(:,21)    ! wptBC
2913                                       ! profile 24 is initial profile (sa)
2914                                       ! profiles 25-29 left empty for initial
2915                                       ! profiles
2916       hom(:,1,30,sr) = sums(:,30)     ! u*2
2917       hom(:,1,31,sr) = sums(:,31)     ! v*2
2918       hom(:,1,32,sr) = sums(:,32)     ! w*2
2919       hom(:,1,33,sr) = sums(:,33)     ! pt*2
2920       hom(:,1,34,sr) = sums(:,34)     ! e*
2921       hom(:,1,35,sr) = sums(:,35)     ! w*2pt*
2922       hom(:,1,36,sr) = sums(:,36)     ! w*pt*2
2923       hom(:,1,37,sr) = sums(:,37)     ! w*e*
2924       hom(:,1,38,sr) = sums(:,38)     ! w*3
2925       hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20_wp )**1.5_wp   ! Sw
2926       hom(:,1,40,sr) = sums(:,40)     ! p
2927       hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
2928       hom(:,1,46,sr) = sums(:,46)     ! w*vpt*
2929       hom(:,1,47,sr) = sums(:,45) + sums(:,46)    ! wvpt
2930       hom(:,1,48,sr) = sums(:,48)     ! w"q" (w"qv")
2931       hom(:,1,49,sr) = sums(:,49)     ! w*q* (w*qv*)
2932       hom(:,1,50,sr) = sums(:,48) + sums(:,49)    ! wq (wqv)
2933       hom(:,1,51,sr) = sums(:,51)     ! w"qv"
2934       hom(:,1,52,sr) = sums(:,52)     ! w*qv*
2935       hom(:,1,53,sr) = sums(:,52) + sums(:,51)    ! wq (wqv)
2936       hom(:,1,54,sr) = sums(:,54)     ! ql
2937       hom(:,1,55,sr) = sums(:,55)     ! w*u*u*/dz
2938       hom(:,1,56,sr) = sums(:,56)     ! w*p*/dz
2939       hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho )/dz
2940       hom(:,1,58,sr) = sums(:,58)     ! u"pt"
2941       hom(:,1,59,sr) = sums(:,59)     ! u*pt*
2942       hom(:,1,60,sr) = sums(:,58) + sums(:,59)    ! upt_t
2943       hom(:,1,61,sr) = sums(:,61)     ! v"pt"
2944       hom(:,1,62,sr) = sums(:,62)     ! v*pt*
2945       hom(:,1,63,sr) = sums(:,61) + sums(:,62)    ! vpt_t
2946       hom(:,1,64,sr) = sums(:,64)     ! rho
2947       hom(:,1,65,sr) = sums(:,65)     ! w"sa"
2948       hom(:,1,66,sr) = sums(:,66)     ! w*sa*
2949       hom(:,1,67,sr) = sums(:,65) + sums(:,66)    ! wsa
2950       hom(:,1,68,sr) = sums(:,68)     ! w*p*
2951       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho
2952       hom(:,1,70,sr) = sums(:,70)     ! q*2
2953       hom(:,1,71,sr) = sums(:,71)     ! prho
2954       hom(:,1,72,sr) = hyp * 1E-4_wp     ! hyp in dbar
2955       hom(:,1,73,sr) = sums(:,73)     ! nr
2956       hom(:,1,74,sr) = sums(:,74)     ! qr
2957       hom(:,1,75,sr) = sums(:,75)     ! qc
2958       hom(:,1,76,sr) = sums(:,76)     ! prr (precipitation rate)
2959                                       ! 77 is initial density profile
2960       hom(:,1,78,sr) = ug             ! ug
2961       hom(:,1,79,sr) = vg             ! vg
2962       hom(:,1,80,sr) = w_subs         ! w_subs
2963
2964       IF ( large_scale_forcing )  THEN
2965          hom(:,1,81,sr) = sums_ls_l(:,0)          ! td_lsa_lpt
2966          hom(:,1,82,sr) = sums_ls_l(:,1)          ! td_lsa_q
2967          IF ( use_subsidence_tendencies )  THEN
2968             hom(:,1,83,sr) = sums_ls_l(:,2)       ! td_sub_lpt
2969             hom(:,1,84,sr) = sums_ls_l(:,3)       ! td_sub_q
2970          ELSE
2971             hom(:,1,83,sr) = sums(:,83)           ! td_sub_lpt
2972             hom(:,1,84,sr) = sums(:,84)           ! td_sub_q
2973          ENDIF
2974          hom(:,1,85,sr) = sums(:,85)              ! td_nud_lpt
2975          hom(:,1,86,sr) = sums(:,86)              ! td_nud_q
2976          hom(:,1,87,sr) = sums(:,87)              ! td_nud_u
2977          hom(:,1,88,sr) = sums(:,88)              ! td_nud_v
2978       END IF
2979
2980       hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
2981                                       ! upstream-parts u_x, u_y, u_z, v_x,
2982                                       ! v_y, usw. (in last but one profile)
2983       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
2984                                       ! u*, w'u', w'v', t* (in last profile)
2985
2986       IF ( max_pr_user > 0 )  THEN    ! user-defined profiles
2987          hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
2988                               sums(:,pr_palm+1:pr_palm+max_pr_user)
2989       ENDIF
2990
2991!
2992!--    Determine the boundary layer height using two different schemes.
2993!--    First scheme: Starting from the Earth's (Ocean's) surface, look for the
2994!--    first relative minimum (maximum) of the total heat flux.
2995!--    The corresponding height is assumed as the boundary layer height, if it
2996!--    is less than 1.5 times the height where the heat flux becomes negative
2997!--    (positive) for the first time.
2998       z_i(1) = 0.0_wp
2999       first = .TRUE.
3000
3001       IF ( ocean )  THEN
3002          DO  k = nzt, nzb+1, -1
3003             IF ( first .AND. hom(k,1,18,sr) < 0.0_wp                          &
3004                .AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp )  THEN
3005                first = .FALSE.
3006                height = zw(k)
3007             ENDIF
3008             IF ( hom(k,1,18,sr) < 0.0_wp  .AND.                               &
3009                  abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND.                        &
3010                  hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
3011                IF ( zw(k) < 1.5_wp * height )  THEN
3012                   z_i(1) = zw(k)
3013                ELSE
3014                   z_i(1) = height
3015                ENDIF
3016                EXIT
3017             ENDIF
3018          ENDDO
3019       ELSE
3020          DO  k = nzb, nzt-1
3021             IF ( first .AND. hom(k,1,18,sr) < 0.0_wp                          &
3022                .AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp )  THEN
3023                first = .FALSE.
3024                height = zw(k)
3025             ENDIF
3026             IF ( hom(k,1,18,sr) < 0.0  .AND. &
3027                  abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND. &
3028                  hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
3029                IF ( zw(k) < 1.5_wp * height )  THEN
3030                   z_i(1) = zw(k)
3031                ELSE
3032                   z_i(1) = height
3033                ENDIF
3034                EXIT
3035             ENDIF
3036          ENDDO
3037       ENDIF
3038
3039!
3040!--    Second scheme: Gradient scheme from Sullivan et al. (1998), modified
3041!--    by Uhlenbrock(2006). The boundary layer height is the height with the
3042!--    maximal local temperature gradient: starting from the second (the last
3043!--    but one) vertical gridpoint, the local gradient must be at least
3044!--    0.2K/100m and greater than the next four gradients.
3045!--    WARNING: The threshold value of 0.2K/100m must be adjusted for the
3046!--             ocean case!
3047       z_i(2) = 0.0_wp
3048       DO  k = nzb+1, nzt+1
3049          dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
3050       ENDDO
3051       dptdz_threshold = 0.2_wp / 100.0_wp
3052
3053       IF ( ocean )  THEN
3054          DO  k = nzt+1, nzb+5, -1
3055             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
3056                  dptdz(k) > dptdz(k-1)  .AND.  dptdz(k) > dptdz(k-2)  .AND.  &
3057                  dptdz(k) > dptdz(k-3)  .AND.  dptdz(k) > dptdz(k-4) )  THEN
3058                z_i(2) = zw(k-1)
3059                EXIT
3060             ENDIF
3061          ENDDO
3062       ELSE
3063          DO  k = nzb+1, nzt-3
3064             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
3065                  dptdz(k) > dptdz(k+1)  .AND.  dptdz(k) > dptdz(k+2)  .AND.  &
3066                  dptdz(k) > dptdz(k+3)  .AND.  dptdz(k) > dptdz(k+4) )  THEN
3067                z_i(2) = zw(k-1)
3068                EXIT
3069             ENDIF
3070          ENDDO
3071       ENDIF
3072
3073       hom(nzb+6,1,pr_palm,sr) = z_i(1)
3074       hom(nzb+7,1,pr_palm,sr) = z_i(2)
3075
3076!
3077!--    Computation of both the characteristic vertical velocity and
3078!--    the characteristic convective boundary layer temperature.
3079!--    The horizontal average at nzb+1 is input for the average temperature.
3080       IF ( hom(nzb,1,18,sr) > 0.0_wp .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8_wp  &
3081           .AND.  z_i(1) /= 0.0_wp )  THEN
3082          hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) *                 &
3083                                       hom(nzb,1,18,sr) *                      &
3084                                       ABS( z_i(1) ) )**0.333333333_wp
3085!--       so far this only works if Prandtl layer is used
3086          hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
3087       ELSE
3088          hom(nzb+8,1,pr_palm,sr)  = 0.0_wp
3089          hom(nzb+11,1,pr_palm,sr) = 0.0_wp
3090       ENDIF
3091
3092!
3093!--    Collect the time series quantities
3094       ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)     ! E
3095       ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)     ! E*
3096       ts_value(3,sr) = dt_3d
3097       ts_value(4,sr) = hom(nzb,1,pr_palm,sr)       ! u*
3098       ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr)     ! th*
3099       ts_value(6,sr) = u_max
3100       ts_value(7,sr) = v_max
3101       ts_value(8,sr) = w_max
3102       ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr)    ! new divergence
3103       ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr)    ! old Divergence
3104       ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr)    ! z_i(1)
3105       ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr)    ! z_i(2)
3106       ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr)    ! w*
3107       ts_value(14,sr) = hom(nzb,1,16,sr)           ! w'pt'   at k=0
3108       ts_value(15,sr) = hom(nzb+1,1,16,sr)         ! w'pt'   at k=1
3109       ts_value(16,sr) = hom(nzb+1,1,18,sr)         ! wpt     at k=1
3110       ts_value(17,sr) = hom(nzb,1,4,sr)            ! pt(0)
3111       ts_value(18,sr) = hom(nzb+1,1,4,sr)          ! pt(zp)
3112       ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr)    ! u'w'    at k=0
3113       ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)    ! v'w'    at k=0
3114       ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
3115
3116       IF ( ts_value(5,sr) /= 0.0_wp )  THEN
3117          ts_value(22,sr) = ts_value(4,sr)**2_wp / &
3118                            ( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
3119       ELSE
3120          ts_value(22,sr) = 10000.0_wp
3121       ENDIF
3122
3123       ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
3124!
3125!--    Calculate additional statistics provided by the user interface
3126       CALL user_statistics( 'time_series', sr, 0 )
3127
3128    ENDDO    ! loop of the subregions
3129
3130    !$acc end data
3131
3132!
3133!-- If required, sum up horizontal averages for subsequent time averaging
3134    IF ( do_sum )  THEN
3135       IF ( average_count_pr == 0 )  hom_sum = 0.0_wp
3136       hom_sum = hom_sum + hom(:,1,:,:)
3137       average_count_pr = average_count_pr + 1
3138       do_sum = .FALSE.
3139    ENDIF
3140
3141!
3142!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
3143!-- This flag is reset after each time step in time_integration.
3144    flow_statistics_called = .TRUE.
3145
3146    CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
3147
3148
3149 END SUBROUTINE flow_statistics
3150#endif
Note: See TracBrowser for help on using the repository browser.