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

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

last commit documented

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