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

Last change on this file since 2101 was 2101, checked in by suehring, 4 years ago

last commit documented

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