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

Last change on this file since 1691 was 1691, checked in by maronga, 6 years ago

various bugfixes and modifications of the atmosphere-land-surface-radiation interaction. Completely re-written routine to calculate surface fluxes (surface_layer_fluxes.f90) that replaces prandtl_fluxes. Minor formatting corrections and renamings

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