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

Last change on this file since 2026 was 2026, checked in by suehring, 5 years ago

Bugfix, enable output of s*2; calculation of domain_averaged perturbation energy; some formatting adjustments

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