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

Last change on this file since 1822 was 1822, checked in by hoffmann, 5 years ago

changes in LPM and bulk cloud microphysics

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