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

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

last commit documented

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