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

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

last commit documented

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