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

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

last commit documented

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