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

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

last commit documented

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