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

Last change on this file since 2073 was 2073, checked in by raasch, 5 years ago

mostly openmp related bugfixes

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