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

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

last commit documented

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