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

Last change on this file since 1658 was 1658, checked in by raasch, 6 years ago

bugfix: temporary reduction variables in openacc brach are now initialized to zero

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