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

Last change on this file since 1353 was 1353, checked in by heinze, 7 years ago

REAL constants provided with KIND-attribute

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