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

Last change on this file since 1965 was 1961, checked in by suehring, 8 years ago

last commit documented

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