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
Line 
1!> @file flow_statistics.f90
2!------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
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.
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!
17! Copyright 1997-2016 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22! Bugfix, enable output of s*2.
23! Change, calculation of domain-averaged perturbation energy.
24! Some formatting adjustments.
25!
26! Former revisions:
27! -----------------
28! $Id: flow_statistics.f90 2026 2016-10-18 10:27:02Z suehring $
29!
30! 2000 2016-08-20 18:09:15Z knoop
31! Forced header and separation lines into 80 columns
32!
33! 1976 2016-07-27 13:28:04Z maronga
34! Removed some unneeded __rrtmg preprocessor directives
35!
36! 1960 2016-07-12 16:34:24Z suehring
37! Separate humidity and passive scalar
38!
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!
43! 1853 2016-04-11 09:00:35Z maronga
44! Adjusted for use with radiation_scheme = constant
45!
46! 1849 2016-04-08 11:33:18Z hoffmann
47! prr moved to arrays_3d
48!
49! 1822 2016-04-07 07:49:42Z hoffmann
50! Output of bulk microphysics simplified.
51!
52! 1815 2016-04-06 13:49:59Z raasch
53! cpp-directives for intel openmp bug removed
54!
55! 1783 2016-03-06 18:36:17Z raasch
56! +module netcdf_interface
57!
58! 1747 2016-02-08 12:25:53Z raasch
59! small bugfixes for accelerator version
60!
61! 1738 2015-12-18 13:56:05Z raasch
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
66!
67! 1709 2015-11-04 14:47:01Z maronga
68! Updated output of Obukhov length
69!
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!
74! 1682 2015-10-07 23:56:08Z knoop
75! Code annotations made doxygen readable
76!
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!
81! 1654 2015-09-17 09:20:17Z raasch
82! FORTRAN bugfix of r1652
83!
84! 1652 2015-09-17 08:12:24Z raasch
85! bugfix in calculation of energy production by turbulent transport of TKE
86!
87! 1593 2015-05-16 13:58:02Z raasch
88! FORTRAN errors removed from openacc branch
89!
90! 1585 2015-04-30 07:05:52Z maronga
91! Added output of timeseries and profiles for RRTMG
92!
93! 1571 2015-03-12 16:12:49Z maronga
94! Bugfix: output of rad_net and rad_sw_in
95!
96! 1567 2015-03-10 17:57:55Z suehring
97! Reverse modifications made for monotonic limiter.
98!
99! 1557 2015-03-05 16:43:04Z suehring
100! Adjustments for monotonic limiter
101!
102! 1555 2015-03-04 17:44:27Z maronga
103! Added output of r_a and r_s.
104!
105! 1551 2015-03-03 14:18:16Z maronga
106! Added suppport for land surface model and radiation model output.
107!
108! 1498 2014-12-03 14:09:51Z suehring
109! Comments added
110!
111! 1482 2014-10-18 12:34:45Z raasch
112! missing ngp_sums_ls added in accelerator version
113!
114! 1450 2014-08-21 07:31:51Z heinze
115! bugfix: calculate fac only for simulated_time >= 0.0
116!
117! 1396 2014-05-06 13:37:41Z raasch
118! bugfix: "copyin" replaced by "update device" in openacc-branch
119!
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!
124! 1382 2014-04-30 12:15:41Z boeske
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
130!
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!
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!
139! 1353 2014-04-08 15:21:23Z heinze
140! REAL constants provided with KIND-attribute
141!
142! 1322 2014-03-20 16:38:49Z raasch
143! REAL constants defined as wp-kind
144!
145! 1320 2014-03-20 08:40:49Z raasch
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
152!
153! 1299 2014-03-06 13:15:21Z heinze
154! Output of large scale vertical velocity w_subs
155!
156! 1257 2013-11-08 15:18:40Z raasch
157! openacc "end parallel" replaced by "end parallel loop"
158!
159! 1241 2013-10-30 11:36:58Z heinze
160! Output of ug and vg
161!
162! 1221 2013-09-10 08:59:13Z raasch
163! ported for openACC in separate #else branch
164!
165! 1179 2013-06-14 05:57:58Z raasch
166! comment for profile 77 added
167!
168! 1115 2013-03-26 18:16:16Z hoffmann
169! ql is calculated by calc_liquid_water_content
170!
171! 1111 2013-03-08 23:54:10Z raasch
172! openACC directive added
173!
174! 1053 2012-11-13 17:11:03Z hoffmann
175! additions for two-moment cloud physics scheme:
176! +nr, qr, qc, prr
177!
178! 1036 2012-10-22 13:43:42Z raasch
179! code put under GPL (PALM 3.9)
180!
181! 1007 2012-09-19 14:30:36Z franke
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
186!
187! 801 2012-01-10 17:30:36Z suehring
188! Calculation of turbulent fluxes in advec_ws is now thread-safe.
189!
190! Revision 1.1  1997/08/11 06:15:17  raasch
191! Initial revision
192!
193!
194! Description:
195! ------------
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.
204!------------------------------------------------------------------------------!
205#if ! defined( __openacc )
206 SUBROUTINE flow_statistics
207 
208
209    USE arrays_3d,                                                             &
210        ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, ol, p, prho, prr, pt, q, qc, ql,&
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
215       
216    USE cloud_parameters,                                                      &
217        ONLY:   l_d_cp, pt_d_t
218       
219    USE control_parameters,                                                    &
220        ONLY:   average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
221                dt_3d, g, humidity, kappa, large_scale_forcing,                &
222                large_scale_subsidence, max_pr_user, message_string, neutral,  &
223                microphysics_seifert, ocean, passive_scalar, simulated_time,   &
224                use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
225                ws_scheme_mom, ws_scheme_sca
226       
227    USE cpulog,                                                                &
228        ONLY:   cpu_log, log_point
229       
230    USE grid_variables,                                                        &
231        ONLY:   ddx, ddy
232       
233    USE indices,                                                               &
234        ONLY:   ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums,      &
235                ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzb_diff_s_inner,        &
236                nzb_s_inner, nzt, nzt_diff
237       
238    USE kinds
239   
240    USE land_surface_model_mod,                                                &
241        ONLY:   ghf_eb, land_surface, m_soil, nzb_soil, nzt_soil,              &
242                qsws_eb, qsws_liq_eb, qsws_soil_eb, qsws_veg_eb, r_a, r_s,     &
243                shf_eb, t_soil
244
245    USE netcdf_interface,                                                      &
246        ONLY:  dots_rad, dots_soil
247
248    USE pegrid
249
250    USE radiation_model_mod,                                                   &
251        ONLY:  radiation, radiation_scheme, rad_net,                           &
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
254
255#if defined ( __rrtmg )
256    USE radiation_model_mod,                                                   &
257        ONLY:  rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
258#endif
259 
260    USE statistics
261
262
263    IMPLICIT NONE
264
265    INTEGER(iwp) ::  i                   !<
266    INTEGER(iwp) ::  j                   !<
267    INTEGER(iwp) ::  k                   !<
268    INTEGER(iwp) ::  k_surface_level     !<
269    INTEGER(iwp) ::  nt                  !<
270    INTEGER(iwp) ::  omp_get_thread_num  !<
271    INTEGER(iwp) ::  sr                  !<
272    INTEGER(iwp) ::  tn                  !<
273   
274    LOGICAL ::  first  !<
275   
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)           !<
290   
291    REAL(wp) ::  dptdz(nzb+1:nzt+1)    !<
292    REAL(wp) ::  sums_ll(nzb:nzt+1,2)  !<
293
294    CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
295
296    !$acc update host( km, kh, e, ol, pt, qs, qsws, shf, ts, u, usws, v, vsws, w )
297
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
302
303       message_string = 'flow_statistics is called two times within one ' // &
304                        'timestep'
305       CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
306
307    ENDIF
308
309!
310!-- Compute statistics for each (sub-)region
311    DO  sr = 0, statistic_regions
312
313!
314!--    Initialize (local) summation array
315       sums_l = 0.0_wp
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
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
325
326!
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.
339!--    Copy the turbulent quantities, evaluated in the advection routines to
340!--    the local array sums_l() for further computations.
341       IF ( ws_scheme_mom .AND. sr == 0 )  THEN
342
343!
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
347             sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
348             sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)
349          ENDIF
350
351          DO  i = 0, threads_per_task-1
352!
353!--          Swap the turbulent quantities evaluated in advec_ws.
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
359             sums_l(:,34,i) = sums_l(:,34,i) + 0.5_wp *                        & 
360                              ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +      &
361                                sums_ws2_ws_l(:,i) )    ! e*
362          ENDDO
363
364       ENDIF
365
366       IF ( ws_scheme_sca .AND. sr == 0 )  THEN
367
368          DO  i = 0, threads_per_task-1
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*
373          ENDDO
374
375       ENDIF
376!
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
381
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
388             DO  k = nzb_s_inner(j,i), nzt+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!
397!--    Horizontally averaged profile of salinity
398       IF ( ocean )  THEN
399          !$OMP DO
400          DO  i = nxl, nxr
401             DO  j =  nys, nyn
402                DO  k = nzb_s_inner(j,i), nzt+1
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!
411!--    Horizontally averaged profiles of virtual potential temperature,
412!--    total water content, specific humidity and liquid water potential
413!--    temperature
414       IF ( humidity )  THEN
415          !$OMP DO
416          DO  i = nxl, nxr
417             DO  j =  nys, nyn
418                DO  k = nzb_s_inner(j,i), nzt+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
430                   DO  k = nzb_s_inner(j,i), nzt+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
448                DO  k = nzb_s_inner(j,i), nzt+1
449                   sums_l(k,117,tn)  = sums_l(k,117,tn) + s(k,j,i) * rmask(j,i,sr)
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)
462             IF ( ocean )  THEN
463                sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)
464             ENDIF
465             IF ( humidity )  THEN
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
474                sums_l(:,117,0) = sums_l(:,117,0) + sums_l(:,117,i)
475             ENDIF
476          ENDDO
477       ENDIF
478
479#if defined( __parallel )
480!
481!--    Compute total sum from local sums
482       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
483       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL,  &
484                           MPI_SUM, comm2d, ierr )
485       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
486       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL,  &
487                           MPI_SUM, comm2d, ierr )
488       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
489       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL,  &
490                           MPI_SUM, comm2d, ierr )
491       IF ( ocean )  THEN
492          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
493          CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb,       &
494                              MPI_REAL, MPI_SUM, comm2d, ierr )
495       ENDIF
496       IF ( humidity ) THEN
497          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
498          CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb,       &
499                              MPI_REAL, MPI_SUM, comm2d, ierr )
500          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
501          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb,       &
502                              MPI_REAL, MPI_SUM, comm2d, ierr )
503          IF ( cloud_physics ) THEN
504             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
505             CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb,    &
506                                 MPI_REAL, MPI_SUM, comm2d, ierr )
507             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
508             CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb,    &
509                                 MPI_REAL, MPI_SUM, comm2d, ierr )
510          ENDIF
511       ENDIF
512
513       IF ( passive_scalar )  THEN
514          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
515          CALL MPI_ALLREDUCE( sums_l(nzb,117,0), sums(nzb,117), nzt+2-nzb,       &
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)
522       IF ( ocean )  sums(:,23) = sums_l(:,23,0)
523       IF ( humidity ) THEN
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
531       IF ( passive_scalar )  sums(:,117) = sums_l(:,117,0)
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.
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)
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
544
545!
546!--    Salinity
547       IF ( ocean )  THEN
548          sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr)
549          hom(:,1,23,sr) = sums(:,23)             ! sa
550       ENDIF
551
552!
553!--    Humidity and cloud parameters
554       IF ( humidity ) THEN
555          sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr)
556          sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
557          hom(:,1,44,sr) = sums(:,44)             ! vpt
558          hom(:,1,41,sr) = sums(:,41)             ! qv (q)
559          IF ( cloud_physics ) THEN
560             sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr)
561             sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr)
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
569       IF ( passive_scalar )  hom(:,1,117,sr) = sums(:,117) /                  &
570            ngp_2dh_s_inner(:,sr)                    ! s
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.
576!--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
577!--    ----  speaking the following k-loop would have to be split up and
578!--          rearranged according to the staggered grid.
579!--          However, this implies no error since staggered velocity components
580!--          are zero at the walls and inside buildings.
581       tn = 0
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 )
585!$     tn = omp_get_thread_num()
586
587       !$OMP DO
588       DO  i = nxl, nxr
589          DO  j =  nys, nyn
590             sums_l_etot = 0.0_wp
591             DO  k = nzb_s_inner(j,i), nzt+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)
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
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
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)
615
616                sums_l_etot  = sums_l_etot + &
617                                        0.5_wp * ( u(k,j,i)**2 + v(k,j,i)**2 + &
618                                        w(k,j,i)**2 ) * rmask(j,i,sr)
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.
625             sums_l(nzb+4,pr_palm,tn) = sums_l(nzb+4,pr_palm,tn) + sums_l_etot
626!
627!--          2D-arrays (being collected in the last column of sums_l)
628             sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +               &
629                                        us(j,i)   * rmask(j,i,sr)
630             sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +             &
631                                        usws(j,i) * rmask(j,i,sr)
632             sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +             &
633                                        vsws(j,i) * rmask(j,i,sr)
634             sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +             &
635                                        ts(j,i)   * rmask(j,i,sr)
636             IF ( humidity )  THEN
637                sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +        &
638                                            qs(j,i)   * rmask(j,i,sr)
639             ENDIF
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
644          ENDDO
645       ENDDO
646
647!
648!--    Computation of statistics when ws-scheme is not used. Else these
649!--    quantities are evaluated in the advection routines.
650       IF ( .NOT. ws_scheme_mom .OR. sr /= 0 .OR. simulated_time == 0.0_wp )   &
651       THEN
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!
666!--                Perturbation energy
667
668                   sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5_wp *                &
669                                  ( ust2 + vst2 + w2 ) * rmask(j,i,sr)
670                ENDDO
671             ENDDO
672          ENDDO
673       ENDIF
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
688
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
695!
696!--    Horizontally averaged profiles of the vertical fluxes
697
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)
704!--          NOTE: for simplicity, nzb_diff_s_inner is used below, although
705!--          ----  strictly speaking the following k-loop would have to be
706!--                split up according to the staggered grid.
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
711!
712!--             Momentum flux w"u"
713                sums_l(k,12,tn) = sums_l(k,12,tn) - 0.25_wp * (                &
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"
721                sums_l(k,14,tn) = sums_l(k,14,tn) - 0.25_wp * (                &
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)                              &
730                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
731                                               * ( pt(k+1,j,i) - pt(k,j,i) )   &
732                                               * ddzu(k+1) * rmask(j,i,sr)
733
734
735!
736!--             Salinity flux w"sa"
737                IF ( ocean )  THEN
738                   sums_l(k,65,tn) = sums_l(k,65,tn)                           &
739                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
740                                               * ( sa(k+1,j,i) - sa(k,j,i) )   &
741                                               * ddzu(k+1) * rmask(j,i,sr)
742                ENDIF
743
744!
745!--             Buoyancy flux, water flux (humidity flux) w"q"
746                IF ( humidity ) THEN
747                   sums_l(k,45,tn) = sums_l(k,45,tn)                           &
748                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
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)                           &
752                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
753                                               * ( q(k+1,j,i) - q(k,j,i) )     &
754                                               * ddzu(k+1) * rmask(j,i,sr)
755
756                   IF ( cloud_physics ) THEN
757                      sums_l(k,51,tn) = sums_l(k,51,tn)                        &
758                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
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
768                   sums_l(k,119,tn) = sums_l(k,119,tn)                         &
769                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
770                                                  * ( s(k+1,j,i) - s(k,j,i) )  &
771                                                  * ddzu(k+1) * rmask(j,i,sr)
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) + &
786                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
787                sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
788                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
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
793                IF ( humidity )  THEN
794                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
795                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
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) *      &
799                                                  qsws(j,i) )
800                   IF ( cloud_droplets )  THEN
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) )
805                   ENDIF
806                   IF ( cloud_physics )  THEN
807!
808!--                   Formula does not work if ql(nzb) /= 0.0
809                      sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                  & 
810                                          qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
811                   ENDIF
812                ENDIF
813                IF ( passive_scalar )  THEN
814                   sums_l(nzb,119,tn) = sums_l(nzb,119,tn) +                     &
815                                        ssws(j,i) * rmask(j,i,sr) ! w"s"
816                ENDIF
817             ENDIF
818
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
825             IF ( land_surface )  THEN
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)
834             ENDIF
835
836             IF ( radiation  .AND.  radiation_scheme /= 'constant' )  THEN
837                sums_l(nzb,101,tn)  = sums_l(nzb,101,tn)  + rad_net(j,i)
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
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)
849                ENDIF
850#endif
851             ENDIF
852!
853!--          Subgridscale fluxes at the top surface
854             IF ( use_top_fluxes )  THEN
855                sums_l(nzt:nzt+1,12,tn) = sums_l(nzt:nzt+1,12,tn) + &
856                                    uswst(j,i) * rmask(j,i,sr)    ! w"u"
857                sums_l(nzt:nzt+1,14,tn) = sums_l(nzt:nzt+1,14,tn) + &
858                                    vswst(j,i) * rmask(j,i,sr)    ! w"v"
859                sums_l(nzt:nzt+1,16,tn) = sums_l(nzt:nzt+1,16,tn) + &
860                                    tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
861                sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + &
862                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
863                sums_l(nzt:nzt+1,61,tn) = sums_l(nzt:nzt+1,61,tn) + &
864                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
865
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
870                IF ( humidity )  THEN
871                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) +                     &
872                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
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) )
877                   IF ( cloud_droplets )  THEN
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) )
882                   ENDIF
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
891                   sums_l(nzt,119,tn) = sums_l(nzt,119,tn) + &
892                                        sswst(j,i) * rmask(j,i,sr) ! w"s"
893                ENDIF
894             ENDIF
895
896!
897!--          Resolved fluxes (can be computed for all horizontal points)
898!--          NOTE: for simplicity, nzb_s_inner is used below, although strictly
899!--          ----  speaking the following k-loop would have to be split up and
900!--                rearranged according to the staggered grid.
901             DO  k = nzb_s_inner(j,i), nzt
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) )
908
909!--             Higher moments
910                sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 *        &
911                                                    rmask(j,i,sr)
912                sums_l(k,36,tn) = sums_l(k,36,tn) + pts**2 * w(k,j,i) *        &
913                                                    rmask(j,i,sr)
914
915!
916!--             Salinity flux and density (density does not belong to here,
917!--             but so far there is no other suitable place to calculate)
918                IF ( ocean )  THEN
919                   IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
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)
924                   ENDIF
925                   sums_l(k,64,tn) = sums_l(k,64,tn) + rho(k,j,i) *            &
926                                                       rmask(j,i,sr)
927                   sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) *           &
928                                                       rmask(j,i,sr)
929                ENDIF
930
931!
932!--             Buoyancy flux, water flux, humidity flux, liquid water
933!--             content, rain drop concentration and rain water content
934                IF ( humidity )  THEN
935                   IF ( cloud_physics .OR. cloud_droplets )  THEN
936                      pts = 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +         &
937                                    vpt(k+1,j,i) - hom(k+1,1,44,sr) )
938                      sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *     &
939                                                          rmask(j,i,sr)
940                      sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * rmask(j,i,sr)
941
942                      IF ( .NOT. cloud_droplets )  THEN
943                         pts = 0.5_wp *                                        &
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) ) -                   &
947                              hom(k+1,1,42,sr) )
948                         sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) *  &
949                                                             rmask(j,i,sr)
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) *    &
956                                                                rmask(j,i,sr)
957                            sums_l(k,74,tn) = sums_l(k,74,tn) + qr(k,j,i) *    &
958                                                                rmask(j,i,sr)
959                         ENDIF
960                      ENDIF
961
962                   ELSE
963                      IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
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) *  &
967                                                             rmask(j,i,sr)
968                      ELSE IF ( ws_scheme_sca .AND. sr == 0 )  THEN
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)
974                      END IF
975                   END IF
976                ENDIF
977!
978!--             Passive scalar flux
979                IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca                &
980                     .OR. sr /= 0 ) )  THEN
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) *      &
984                                                       rmask(j,i,sr)
985                ENDIF
986
987!
988!--             Energy flux w*e*
989!--             has to be adjusted
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 ) &
992                                             * rmask(j,i,sr)
993             ENDDO
994          ENDDO
995       ENDDO
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:
1000       IF ( .NOT. ws_scheme_mom .OR. sr /= 0  )  THEN
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
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) )
1009!
1010!--               Momentum flux w*u*
1011                  sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5_wp *                 &
1012                                                    ( w(k,j,i-1) + w(k,j,i) )  &
1013                                                    * ust * rmask(j,i,sr)
1014!
1015!--               Momentum flux w*v*
1016                  sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5_wp *                 &
1017                                                    ( w(k,j-1,i) + w(k,j,i) )  &
1018                                                    * vst * rmask(j,i,sr)
1019               ENDDO
1020            ENDDO
1021         ENDDO
1022
1023       ENDIF
1024       IF ( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
1025         !$OMP DO
1026         DO  i = nxl, nxr
1027            DO  j = nys, nyn
1028               DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
1029!
1030!--               Vertical heat flux
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) )                   &
1034                           * w(k,j,i) * rmask(j,i,sr)
1035                  IF ( humidity )  THEN
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)
1040                  ENDIF
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) )
1044                     sums_l(k,116,tn) = sums_l(k,116,tn) + pts * w(k,j,i) *     &
1045                                        rmask(j,i,sr)
1046                  ENDIF
1047               ENDDO
1048            ENDDO
1049         ENDDO
1050
1051       ENDIF
1052
1053!
1054!--    Density at top follows Neumann condition
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
1059
1060!
1061!--    Divergence of vertical flux of resolved scale energy and pressure
1062!--    fluctuations as well as flux of pressure fluctuation itself (68).
1063!--    First calculate the products, then the divergence.
1064!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
1065       IF ( hom(nzb+1,2,55,0) /= 0.0_wp  .OR.  hom(nzb+1,2,68,0) /= 0.0_wp )   &
1066       THEN
1067          sums_ll = 0.0_wp  ! local array
1068
1069          !$OMP DO
1070          DO  i = nxl, nxr
1071             DO  j = nys, nyn
1072                DO  k = nzb_s_inner(j,i)+1, nzt
1073
1074                   sums_ll(k,1) = sums_ll(k,1) + 0.5_wp * w(k,j,i) * (         &
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) )  &
1078                            - 0.5_wp * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) ) )**2&
1079                + w(k,j,i)**2                                        )
1080
1081                   sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i)             &
1082                                               * ( p(k,j,i) + p(k+1,j,i) )
1083
1084                ENDDO
1085             ENDDO
1086          ENDDO
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
1091
1092          DO  k = nzb+1, nzt
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)
1095             sums_l(k,68,tn) = sums_ll(k,2)
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)
1099          sums_l(nzb,68,tn) = 0.0_wp    ! because w* = 0 at nzb
1100
1101       ENDIF
1102
1103!
1104!--    Divergence of vertical flux of SGS TKE and the flux itself (69)
1105       IF ( hom(nzb+1,2,57,0) /= 0.0_wp  .OR.  hom(nzb+1,2,69,0) /= 0.0_wp )   &
1106       THEN
1107          !$OMP DO
1108          DO  i = nxl, nxr
1109             DO  j = nys, nyn
1110                DO  k = nzb_s_inner(j,i)+1, nzt
1111
1112                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5_wp * (              &
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)   &
1115                                                                ) * ddzw(k)
1116
1117                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5_wp * (              &
1118                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
1119                                                                )
1120
1121                ENDDO
1122             ENDDO
1123          ENDDO
1124          sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
1125          sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
1126
1127       ENDIF
1128
1129!
1130!--    Horizontal heat fluxes (subgrid, resolved, total).
1131!--    Do it only, if profiles shall be plotted.
1132       IF ( hom(nzb+1,2,58,0) /= 0.0_wp ) THEN
1133
1134          !$OMP DO
1135          DO  i = nxl, nxr
1136             DO  j = nys, nyn
1137                DO  k = nzb_s_inner(j,i)+1, nzt
1138!
1139!--                Subgrid horizontal heat fluxes u"pt", v"pt"
1140                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5_wp *                &
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)
1144                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
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) ) &
1152                                    * 0.5_wp * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
1153                                                 pt(k,j,i)   - hom(k,1,4,sr) )
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) )
1156                   sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
1157                                                  ( v(k,j,i) - hom(k,1,2,sr) ) &
1158                                    * 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
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)
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
1171
1172       ENDIF
1173
1174!
1175!--    Collect current large scale advection and subsidence tendencies for
1176!--    data output
1177       IF ( large_scale_forcing  .AND.  ( simulated_time > 0.0_wp ) )  THEN
1178!
1179!--       Interpolation in time of LSF_DATA
1180          nt = 1
1181          DO WHILE ( simulated_time - dt_3d > time_vert(nt) )
1182             nt = nt + 1
1183          ENDDO
1184          IF ( simulated_time - dt_3d /= time_vert(nt) )  THEN
1185            nt = nt - 1
1186          ENDIF
1187
1188          fac = ( simulated_time - dt_3d - time_vert(nt) )                     &
1189                / ( time_vert(nt+1)-time_vert(nt) )
1190
1191
1192          DO  k = nzb, nzt
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) )
1197          ENDDO
1198
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
1202          IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
1203
1204             DO  k = nzb, nzt
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) )
1209             ENDDO
1210
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
1214          ENDIF
1215
1216       ENDIF
1217
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
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)
1228                ENDDO
1229             ENDDO
1230          ENDDO
1231       ENDIF
1232       
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
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)
1252                   sums_l(k,109,tn)  = sums_l(k,109,tn)  + rad_sw_hr(k,j,i)    &
1253                                       * rmask(j,i,sr)
1254                ENDDO
1255             ENDDO
1256          ENDDO
1257       ENDIF
1258!
1259!--    Calculate the user-defined profiles
1260       CALL user_statistics( 'profiles', sr, tn )
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)
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
1276          ENDDO
1277       ENDIF
1278
1279#if defined( __parallel )
1280
1281!
1282!--    Compute total sum from local sums
1283       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1284       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL,   &
1285                           MPI_SUM, comm2d, ierr )
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
1290#else
1291       sums = sums_l(:,:,0)
1292       IF ( large_scale_forcing )  THEN
1293          sums(:,81:88) = sums_ls_l
1294       ENDIF
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.
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
1303!--    Profiles:
1304       DO  k = nzb, nzt+1
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)
1313          sums(k,116)           = sums(k,116)           / ngp_2dh(sr)
1314          sums(k,119)           = sums(k,119)           / ngp_2dh(sr)
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)
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)
1325          ENDIF
1326       ENDDO
1327
1328!--    u* and so on
1329!--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
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)
1332       sums(nzb:nzb+3,pr_palm)    = sums(nzb:nzb+3,pr_palm)    / &
1333                                    ngp_2dh(sr)
1334       sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)       / &    ! qs
1335                                    ngp_2dh(sr)
1336       sums(nzb+13,pr_palm)       = sums(nzb+13,pr_palm)       / &    ! ss
1337                                    ngp_2dh(sr)
1338!--    eges, e*
1339       sums(nzb+4:nzb+5,pr_palm)  = sums(nzb+4:nzb+5,pr_palm)  / &
1340                                    ngp_3d(sr)
1341!--    Old and new divergence
1342       sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
1343                                    ngp_3d_inner(sr)
1344
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) / &
1350                                    ngp_2dh_s_inner(k,sr)
1351          ENDDO
1352       ENDIF
1353
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
1373                                       ! profile 24 is initial profile (sa)
1374                                       ! profiles 25-29 left empty for initial
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
1385       hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20_wp )**1.5_wp   ! Sw
1386       hom(:,1,40,sr) = sums(:,40)     ! p
1387       hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
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
1399       hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho )/dz
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
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
1410       hom(:,1,68,sr) = sums(:,68)     ! w*p*
1411       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho
1412       hom(:,1,70,sr) = sums(:,70)     ! q*2
1413       hom(:,1,71,sr) = sums(:,71)     ! prho
1414       hom(:,1,72,sr) = hyp * 1E-4_wp  ! hyp in dbar
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)
1419                                       ! 77 is initial density profile
1420       hom(:,1,78,sr) = ug             ! ug
1421       hom(:,1,79,sr) = vg             ! vg
1422       hom(:,1,80,sr) = w_subs         ! w_subs
1423
1424       IF ( large_scale_forcing )  THEN
1425          hom(:,1,81,sr) = sums_ls_l(:,0)          ! td_lsa_lpt
1426          hom(:,1,82,sr) = sums_ls_l(:,1)          ! td_lsa_q
1427          IF ( use_subsidence_tendencies )  THEN
1428             hom(:,1,83,sr) = sums_ls_l(:,2)       ! td_sub_lpt
1429             hom(:,1,84,sr) = sums_ls_l(:,3)       ! td_sub_q
1430          ELSE
1431             hom(:,1,83,sr) = sums(:,83)           ! td_sub_lpt
1432             hom(:,1,84,sr) = sums(:,84)           ! td_sub_q
1433          ENDIF
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
1438       ENDIF
1439
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
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
1454       ENDIF
1455
1456       IF ( radiation )  THEN
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
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
1473          ENDIF
1474       ENDIF
1475
1476       hom(:,1,114,sr) = sums(:,114)            !: L
1477
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
1482          hom(:,1,118,sr) = sums(:,118)     ! s*2
1483       ENDIF
1484
1485       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
1486                                       ! u*, w'u', w'v', t* (in last profile)
1487
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
1493!
1494!--    Determine the boundary layer height using two different schemes.
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.
1500       z_i(1) = 0.0_wp
1501       first = .TRUE.
1502
1503       IF ( ocean )  THEN
1504          DO  k = nzt, nzb+1, -1
1505             IF ( first  .AND.  hom(k,1,18,sr) < -1.0E-8_wp )  THEN
1506                first = .FALSE.
1507                height = zw(k)
1508             ENDIF
1509             IF ( hom(k,1,18,sr) < -1.0E-8_wp  .AND.                           &
1510                  hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
1511                IF ( zw(k) < 1.5_wp * height )  THEN
1512                   z_i(1) = zw(k)
1513                ELSE
1514                   z_i(1) = height
1515                ENDIF
1516                EXIT
1517             ENDIF
1518          ENDDO
1519       ELSE
1520          DO  k = nzb, nzt-1
1521             IF ( first  .AND.  hom(k,1,18,sr) < -1.0E-8_wp )  THEN
1522                first = .FALSE.
1523                height = zw(k)
1524             ENDIF
1525             IF ( hom(k,1,18,sr) < -1.0E-8_wp  .AND.                           &
1526                  hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
1527                IF ( zw(k) < 1.5_wp * height )  THEN
1528                   z_i(1) = zw(k)
1529                ELSE
1530                   z_i(1) = height
1531                ENDIF
1532                EXIT
1533             ENDIF
1534          ENDDO
1535       ENDIF
1536
1537!
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!
1545       z_i(2) = 0.0_wp
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
1549       dptdz_threshold = 0.2_wp / 100.0_wp
1550
1551       IF ( ocean )  THEN
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)
1557                EXIT
1558             ENDIF
1559          ENDDO
1560       ELSE
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)
1566                EXIT
1567             ENDIF
1568          ENDDO
1569       ENDIF
1570
1571       hom(nzb+6,1,pr_palm,sr) = z_i(1)
1572       hom(nzb+7,1,pr_palm,sr) = z_i(2)
1573
1574!
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!
1584!--    Computation of both the characteristic vertical velocity and
1585!--    the characteristic convective boundary layer temperature.
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
1595       ELSE
1596          hom(nzb+8,1,pr_palm,sr)  = 0.0_wp
1597       ENDIF
1598
1599!
1600!--    Collect the time series quantities
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*
1603       ts_value(3,sr) = dt_3d
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*
1606       ts_value(6,sr) = u_max
1607       ts_value(7,sr) = v_max
1608       ts_value(8,sr) = w_max
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*
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)
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
1621       ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
1622
1623       IF ( .NOT. neutral )  THEN
1624          ts_value(22,sr) = hom(nzb,1,114,sr)          ! L
1625       ELSE
1626          ts_value(22,sr) = 1.0E10_wp
1627       ENDIF
1628
1629       ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
1630
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
1636!
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
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
1647       ENDIF
1648!
1649!--    Collect radiation model timeseries
1650       IF ( radiation )  THEN
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
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
1656
1657          IF ( radiation_scheme == 'rrtmg' )  THEN
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
1662          ENDIF
1663
1664       ENDIF
1665
1666!
1667!--    Calculate additional statistics provided by the user interface
1668       CALL user_statistics( 'time_series', sr, 0 )
1669
1670    ENDDO    ! loop of the subregions
1671
1672!
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
1676       IF ( average_count_pr == 0 )  hom_sum = 0.0_wp
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
1691
1692
1693#else
1694
1695
1696!------------------------------------------------------------------------------!
1697! Description:
1698! ------------
1699!> flow statistics - accelerator version
1700!------------------------------------------------------------------------------!
1701 SUBROUTINE flow_statistics
1702
1703    USE arrays_3d,                                                             &
1704        ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, p, prho, pt, q, qc, ql, qr, qs, &
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
1709                 
1710       
1711    USE cloud_parameters,                                                      &
1712        ONLY:  l_d_cp, prr, pt_d_t
1713       
1714    USE control_parameters,                                                    &
1715        ONLY :  average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
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
1721       
1722    USE cpulog,                                                                &
1723        ONLY:  cpu_log, log_point
1724       
1725    USE grid_variables,                                                        &
1726        ONLY:  ddx, ddy
1727       
1728    USE indices,                                                               &
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
1732       
1733    USE kinds
1734   
1735    USE land_surface_model_mod,                                                &
1736        ONLY:   ghf_eb, land_surface, m_soil, nzb_soil, nzt_soil,              &
1737                qsws_eb, qsws_liq_eb, qsws_soil_eb, qsws_veg_eb, r_a, r_s,     &
1738                shf_eb, t_soil
1739
1740    USE netcdf_interface,                                                      &
1741        ONLY:  dots_rad, dots_soil
1742
1743    USE pegrid
1744   
1745    USE radiation_model_mod,                                                   &
1746        ONLY:  radiation, radiation_scheme, rad_net,                 &
1747               rad_lw_in, rad_lw_out, rad_sw_in, rad_sw_out
1748
1749#if defined ( __rrtmg )
1750    USE radiation_model_mod,                                                   &
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
1753#endif
1754
1755    USE statistics
1756
1757    IMPLICIT NONE
1758
1759    INTEGER(iwp) ::  i                   !<
1760    INTEGER(iwp) ::  j                   !<
1761    INTEGER(iwp) ::  k                   !<
1762    INTEGER(iwp) ::  k_surface_level     !<
1763    INTEGER(iwp) ::  nt                  !<
1764    INTEGER(iwp) ::  omp_get_thread_num  !<
1765    INTEGER(iwp) ::  sr                  !<
1766    INTEGER(iwp) ::  tn                  !<
1767   
1768    LOGICAL ::  first  !<
1769   
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)           !<
1791
1792    REAL(wp) ::  dptdz(nzb+1:nzt+1)    !<
1793    REAL(wp) ::  sums_ll(nzb:nzt+1,2)  !<
1794
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
1808    !$acc data create( sums, sums_l )
1809    !$acc update device( hom )
1810
1811!
1812!-- Compute statistics for each (sub-)region
1813    DO  sr = 0, statistic_regions
1814
1815!
1816!--    Initialize (local) summation array
1817       sums_l = 0.0_wp
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!
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.
1841!--    Copy the turbulent quantities, evaluated in the advection routines to
1842!--    the local array sums_l() for further computations.
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
1861             sums_l(:,34,i) = sums_l(:,34,i) + 0.5_wp *                        &
1862                              ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +      &
1863                                sums_ws2_ws_l(:,i) )    ! e*
1864             DO  k = nzb, nzt
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)     )
1869             ENDDO
1870          ENDDO
1871
1872       ENDIF
1873
1874       IF ( ws_scheme_sca .AND. sr == 0 )  THEN
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*
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*
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
1898          s1 = 0
1899          s2 = 0
1900          s3 = 0
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
1915       !$acc end parallel loop
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
1923             s1 = 0
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
1932          !$acc end parallel loop
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
1944             s1 = 0
1945             s2 = 0
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
1956          !$acc end parallel loop
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
1962                s1 = 0
1963                s2 = 0
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
1976             !$acc end parallel loop
1977          ENDIF
1978       ENDIF
1979
1980!
1981!--    Horizontally averaged profiles of passive scalar
1982       IF ( passive_scalar )  THEN
1983          !$OMP DO
1984          !$acc parallel loop gang present( s, rflags_invers, rmask, sums_l ) create( s1 )
1985          DO  k = nzb, nzt+1
1986             s1 = 0
1987             !$acc loop vector collapse( 2 ) reduction( +: s1 )
1988             DO  i = nxl, nxr
1989                DO  j =  nys, nyn
1990                   s1 = s1 + s(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1991                ENDDO
1992             ENDDO
1993             sums_l(k,117,tn) = s1
1994          ENDDO
1995          !$acc end parallel loop
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 )
2027                sums_l(:,117,0) = sums_l(:,117,0) + sums_l(:,117,i)
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 )
2038       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL,  &
2039                           MPI_SUM, comm2d, ierr )
2040       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2041       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL,  &
2042                           MPI_SUM, comm2d, ierr )
2043       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2044       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL,  &
2045                           MPI_SUM, comm2d, ierr )
2046       IF ( ocean )  THEN
2047          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2048          CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb,       &
2049                              MPI_REAL, MPI_SUM, comm2d, ierr )
2050       ENDIF
2051       IF ( humidity ) THEN
2052          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2053          CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb,       &
2054                              MPI_REAL, MPI_SUM, comm2d, ierr )
2055          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2056          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb,       &
2057                              MPI_REAL, MPI_SUM, comm2d, ierr )
2058          IF ( cloud_physics ) THEN
2059             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2060             CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb,    &
2061                                 MPI_REAL, MPI_SUM, comm2d, ierr )
2062             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2063             CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb,    &
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 )
2070          CALL MPI_ALLREDUCE( sums_l(nzb,117,0), sums(nzb,117), nzt+2-nzb,     &
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 )
2099          sums(:,117) = sums_l(:,117,0)
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 )
2148          sums(:,117)     = sums(:,117) / ngp_2dh_s_inner(:,sr)
2149          hom(:,1,117,sr) = sums(:,117)                ! s
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
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 )
2166!$     tn = omp_get_thread_num()
2167
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
2171          s1 = 0
2172          s2 = 0
2173          s3 = 0
2174          s4 = 0
2175          s5 = 0
2176          s6 = 0
2177          s7 = 0
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
2204       !$acc end parallel loop
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
2210             s1 = 0
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
2220          !$acc end parallel loop
2221       ENDIF
2222
2223!
2224!--    Total and perturbation energy for the total domain (being
2225!--    collected in the last column of sums_l).
2226       s1 = 0
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
2232                s1 = s1 + 0.5_wp *                                             & 
2233                          ( u(k,j,i)**2 + v(k,j,i)**2 + w(k,j,i)**2 ) *        &
2234                          rmask(j,i,sr) * rflags_invers(j,i,k+1)
2235             ENDDO
2236          ENDDO
2237       ENDDO
2238       !$acc end parallel loop
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 )
2245       s1 = 0
2246       s2 = 0
2247       s3 = 0
2248       s4 = 0
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 )
2268          s1 = 0
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
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
2292!
2293!--    Computation of statistics when ws-scheme is not used. Else these
2294!--    quantities are evaluated in the advection routines.
2295       IF ( .NOT. ws_scheme_mom .OR. sr /= 0 .OR. simulated_time == 0.0_wp )   &
2296       THEN
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
2301             s1 = 0
2302             s2 = 0
2303             s3 = 0
2304             s4 = 0
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
2317                   s4 = s4 + 0.5_wp * ( ust2 + vst2 + w2 ) * rmask(j,i,sr) *   &
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
2326          !$acc end parallel loop
2327!
2328!--       Total perturbation TKE
2329          !$OMP DO
2330          !$acc parallel present( sums_l ) create( s1 )
2331          s1 = 0
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
2359          s1 = 0
2360          s2 = 0
2361          s3 = 0
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"
2368                s1 = s1 - 0.25_wp * (                                          &
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"
2377                s2 = s2 - 0.25_wp * (                                          &
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"
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)
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
2396       !$acc end parallel loop
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
2403             s1 = 0
2404             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2405             DO  i = nxl, nxr
2406                DO  j = nys, nyn
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)
2411                ENDDO
2412             ENDDO
2413             sums_l(k,65,tn) = s1
2414          ENDDO
2415          !$acc end parallel loop
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
2424             s1 = 0
2425             s2 = 0
2426             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2427             DO  i = nxl, nxr
2428                DO  j = nys, nyn
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) )            &
2431                                    * ddzu(k+1) * rmask(j,i,sr)                &
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)
2437                ENDDO
2438             ENDDO
2439             sums_l(k,45,tn) = s1
2440             sums_l(k,48,tn) = s2
2441          ENDDO
2442          !$acc end parallel loop
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
2448                s1 = 0
2449                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2450                DO  i = nxl, nxr
2451                   DO  j = nys, nyn
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)
2457                   ENDDO
2458                ENDDO
2459                sums_l(k,51,tn) = s1
2460             ENDDO
2461             !$acc end parallel loop
2462
2463          ENDIF
2464
2465       ENDIF
2466!
2467!--    Passive scalar flux
2468       IF ( passive_scalar )  THEN
2469
2470          !$acc parallel loop gang present( ddzu, kh, s, rflags_invers, rmask, sums_l ) create( s1 )
2471          DO  k = nzb, nzt_diff
2472             s1 = 0
2473             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2474             DO  i = nxl, nxr
2475                DO  j = nys, nyn
2476                   s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
2477                                    * ( s(k+1,j,i) - s(k,j,i) )                &
2478                                    * ddzu(k+1) * rmask(j,i,sr)                &
2479                                    * rflags_invers(j,i,k+1)
2480                ENDDO
2481             ENDDO
2482             sums_l(k,119,tn) = s1
2483          ENDDO
2484          !$acc end parallel loop
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 )
2492          s1 = 0
2493          s2 = 0
2494          s3 = 0
2495          s4 = 0
2496          s5 = 0
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"
2505                s4 = s4 + 0.0_wp * rmask(j,i,sr)        ! u"pt"
2506                s5 = s5 + 0.0_wp * rmask(j,i,sr)        ! v"pt"
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 )
2520             s1 = 0
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 )
2536             s1 = 0
2537             s2 = 0
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")
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) )
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 )
2554                s1 = 0
2555                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2556                DO  i = nxl, nxr
2557                   DO  j =  nys, nyn
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) )
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 )
2572                s1 = 0
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
2591             !$acc parallel present( ssws, rmask, sums_l ) create( s1 )
2592             s1 = 0
2593             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2594             DO  i = nxl, nxr
2595                DO  j =  nys, nyn
2596                   s1 = s1 + ssws(j,i) * rmask(j,i,sr)  ! w"s"
2597                ENDDO
2598             ENDDO
2599             sums_l(nzb,119,tn) = s1
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 )
2612          s1 = 0
2613          s2 = 0
2614          s3 = 0
2615          s4 = 0
2616          s5 = 0
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"
2623                s4 = s4 + 0.0_wp * rmask(j,i,sr)        ! u"pt"
2624                s5 = s5 + 0.0_wp * rmask(j,i,sr)        ! v"pt"
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 )
2638             s1 = 0
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 )
2654             s1 = 0
2655             s2 = 0
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")
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) )
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 )
2672                s1 = 0
2673                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2674                DO  i = nxl, nxr
2675                   DO  j =  nys, nyn
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) )
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 )
2691                s1 = 0
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
2710             !$acc parallel present( sswst, rmask, sums_l ) create( s1 )
2711             s1 = 0
2712             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2713             DO  i = nxl, nxr
2714                DO  j =  nys, nyn
2715                   s1 = s1 + sswst(j,i) * rmask(j,i,sr) ! w"s"
2716                ENDDO
2717             ENDDO
2718             sums_l(nzt,119,tn) = s1
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
2732          s1 = 0
2733          s2 = 0
2734          s3 = 0
2735          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
2736          DO  i = nxl, nxr
2737             DO  j = nys, nyn
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) )
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?)
2750                s3 = s3 + w(k,j,i) * 0.5_wp * ( ust**2 + vst**2 + w(k,j,i)**2 )&
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
2758       !$acc end parallel loop
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
2765          IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
2766
2767             !$acc parallel loop gang present( hom, rflags_invers, rmask, sa, sums_l, w ) create( s1 )
2768             DO  k = nzb, nzt_diff
2769                s1 = 0
2770                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2771                DO  i = nxl, nxr
2772                   DO  j = nys, nyn
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)
2777                   ENDDO
2778                ENDDO
2779                sums_l(k,66,tn) = s1
2780             ENDDO
2781             !$acc end parallel loop
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
2787             s1 = 0
2788             s2 = 0
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
2799          !$acc end parallel loop
2800
2801       ENDIF
2802
2803!
2804!--    Buoyancy flux, water flux, humidity flux, liquid water
2805!--    content, rain drop concentration and rain water content
2806       IF ( humidity )  THEN
2807
2808          IF ( cloud_physics  .OR.  cloud_droplets )  THEN
2809
2810             !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, vpt, w ) create( s1 )
2811             DO  k = nzb, nzt_diff
2812                s1 = 0
2813                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2814                DO  i = nxl, nxr
2815                   DO  j = nys, nyn
2816                      s1 = s1 + 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +     &
2817                                           vpt(k+1,j,i) - hom(k+1,1,44,sr) ) * &
2818                                         w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2819                   ENDDO
2820                ENDDO
2821                sums_l(k,46,tn) = s1
2822             ENDDO
2823             !$acc end parallel loop
2824
2825             IF ( .NOT. cloud_droplets )  THEN
2826
2827                !$acc parallel loop gang present( hom, q, ql, rflags_invers, rmask, sums_l, w ) create( s1 )
2828                DO  k = nzb, nzt_diff
2829                   s1 = 0
2830                   !$acc loop vector collapse( 2 ) reduction( +: s1 )
2831                   DO  i = nxl, nxr
2832                      DO  j = nys, nyn
2833                         s1 = s1 + 0.5_wp * ( ( q(k,j,i)   - ql(k,j,i)   ) - hom(k,1,42,sr) +   &
2834                                              ( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) ) &
2835                                          * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2836                      ENDDO
2837                   ENDDO
2838                   sums_l(k,52,tn) = s1
2839                ENDDO
2840                !$acc end parallel loop
2841
2842                IF ( microphysics_seifert )  THEN
2843
2844                   !$acc parallel loop gang present( qc, ql, rflags_invers, rmask, sums_l ) create( s1, s2 )
2845                   DO  k = nzb, nzt_diff
2846                      s1 = 0
2847                      !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2848                      DO  i = nxl, nxr
2849                         DO  j = nys, nyn
2850                            s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2851                            s2 = s2 + qc(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2852                         ENDDO
2853                      ENDDO
2854                      sums_l(k,54,tn) = s1
2855                      sums_l(k,75,tn) = s2
2856                   ENDDO
2857                   !$acc end parallel loop
2858
2859                   !$acc parallel loop gang present( nr, qr, prr, rflags_invers, rmask, sums_l ) create( s1, s2, s3 )
2860                   DO  k = nzb, nzt_diff
2861                      s1 = 0
2862                      s2 = 0
2863                      s3 = 0
2864                      !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
2865                      DO  i = nxl, nxr
2866                         DO  j = nys, nyn
2867                            s1 = s1 + nr(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2868                            s2 = s2 + qr(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2869                            s3 = s3 + prr(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2870                         ENDDO
2871                      ENDDO
2872                      sums_l(k,73,tn) = s1
2873                      sums_l(k,74,tn) = s2
2874                      sums_l(k,76,tn) = s3
2875                   ENDDO
2876                   !$acc end parallel loop
2877
2878                ELSE
2879
2880                   !$acc parallel loop gang present( ql, rflags_invers, rmask, sums_l ) create( s1 )
2881                   DO  k = nzb, nzt_diff
2882                      s1 = 0
2883                      !$acc loop vector collapse( 2 ) reduction( +: s1 )
2884                      DO  i = nxl, nxr
2885                         DO  j = nys, nyn
2886                            s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2887                         ENDDO
2888                      ENDDO
2889                      sums_l(k,54,tn) = s1
2890                   ENDDO
2891                   !$acc end parallel loop
2892
2893                ENDIF
2894
2895             ELSE
2896
2897                !$acc parallel loop gang present( ql, rflags_invers, rmask, sums_l ) create( s1 )
2898                DO  k = nzb, nzt_diff
2899                   s1 = 0
2900                   !$acc loop vector collapse( 2 ) reduction( +: s1 )
2901                   DO  i = nxl, nxr
2902                      DO  j = nys, nyn
2903                         s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2904                      ENDDO
2905                   ENDDO
2906                   sums_l(k,54,tn) = s1
2907                ENDDO
2908                !$acc end parallel loop
2909
2910             ENDIF
2911
2912          ELSE
2913
2914             IF( .NOT. ws_scheme_sca  .OR.  sr /= 0 )  THEN
2915
2916                !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, vpt, w ) create( s1 )
2917                DO  k = nzb, nzt_diff
2918                   s1 = 0
2919                   !$acc loop vector collapse( 2 ) reduction( +: s1 )
2920                   DO  i = nxl, nxr
2921                      DO  j = nys, nyn
2922                         s1 = s1 + 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +   &
2923                                              vpt(k+1,j,i) - hom(k+1,1,44,sr) ) &
2924                                          * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2925                      ENDDO
2926                   ENDDO
2927                   sums_l(k,46,tn) = s1
2928                ENDDO
2929                !$acc end parallel loop
2930
2931             ELSEIF ( ws_scheme_sca  .AND.  sr == 0 )  THEN
2932
2933                !$acc parallel loop present( hom, sums_l )
2934                DO  k = nzb, nzt_diff
2935                   sums_l(k,46,tn) = ( 1.0_wp + 0.61_wp * hom(k,1,41,sr) ) * sums_l(k,17,tn) + &
2936                                                0.61_wp * hom(k,1,4,sr) * sums_l(k,49,tn)
2937                ENDDO
2938                !$acc end parallel loop
2939
2940             ENDIF
2941
2942          ENDIF
2943
2944       ENDIF
2945!
2946!--    Passive scalar flux
2947       IF ( passive_scalar  .AND.  ( .NOT. ws_scheme_sca  .OR.  sr /= 0 ) )  THEN
2948
2949          !$acc parallel loop gang present( hom, s, rflags_invers, rmask, sums_l, w ) create( s1 )
2950          DO  k = nzb, nzt_diff
2951             s1 = 0
2952             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2953             DO  i = nxl, nxr
2954                DO  j = nys, nyn
2955                   s1 = s1 + 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +          &
2956                                        s(k+1,j,i) - hom(k+1,1,117,sr) )        &
2957                                    * w(k,j,i) * rmask(j,i,sr)                 &
2958                                    * rflags_invers(j,i,k+1)
2959                ENDDO
2960             ENDDO
2961             sums_l(k,49,tn) = s1
2962          ENDDO
2963          !$acc end parallel loop
2964
2965       ENDIF
2966
2967!
2968!--    For speed optimization fluxes which have been computed in part directly
2969!--    inside the WS advection routines are treated seperatly
2970!--    Momentum fluxes first:
2971       IF ( .NOT. ws_scheme_mom  .OR.  sr /= 0  )  THEN
2972
2973          !$OMP DO
2974          !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, u, v, w ) create( s1, s2 )
2975          DO  k = nzb, nzt_diff
2976             s1 = 0
2977             s2 = 0
2978             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2979             DO  i = nxl, nxr
2980                DO  j = nys, nyn
2981                   ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +               &
2982                                    u(k+1,j,i) - hom(k+1,1,1,sr) )
2983                   vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +               &
2984                                    v(k+1,j,i) - hom(k+1,1,2,sr) )
2985!
2986!--                Momentum flux w*u*
2987                   s1 = s1 + 0.5_wp * ( w(k,j,i-1) + w(k,j,i) )                &
2988                                    * ust * rmask(j,i,sr)                      &
2989                                    * rflags_invers(j,i,k+1)
2990!
2991!--                Momentum flux w*v*
2992                   s2 = s2 + 0.5_wp * ( w(k,j-1,i) + w(k,j,i) )                &
2993                                    * vst * rmask(j,i,sr)                      & 
2994                                    * rflags_invers(j,i,k+1)
2995                ENDDO
2996             ENDDO
2997             sums_l(k,13,tn) = s1
2998             sums_l(k,15,tn) = s1
2999          ENDDO
3000          !$acc end parallel loop
3001
3002       ENDIF
3003
3004       IF ( .NOT. ws_scheme_sca  .OR.  sr /= 0 )  THEN
3005
3006          !$OMP DO
3007          !$acc parallel loop gang present( hom, pt, rflags_invers, rmask, sums_l, w ) create( s1 )
3008          DO  k = nzb, nzt_diff
3009             s1 = 0
3010             !$acc loop vector collapse( 2 ) reduction( +: s1 )
3011             DO  i = nxl, nxr
3012                DO  j = nys, nyn
3013!
3014!--                Vertical heat flux
3015                   s1 = s1 + 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) +          &
3016                                        pt(k+1,j,i) - hom(k+1,1,4,sr) )        &
3017                                    * w(k,j,i) * rmask(j,i,sr)                 & 
3018                                    * rflags_invers(j,i,k+1)
3019                ENDDO
3020             ENDDO
3021             sums_l(k,17,tn) = s1
3022          ENDDO
3023          !$acc end parallel loop
3024
3025          IF ( humidity )  THEN
3026
3027             !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l, w ) create( s1 )
3028             DO  k = nzb, nzt_diff
3029                s1 = 0
3030                !$acc loop vector collapse( 2 ) reduction( +: s1 )
3031                DO  i = nxl, nxr
3032                   DO  j = nys, nyn
3033                      s1 = s1 + 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +       &
3034                                           q(k+1,j,i) - hom(k+1,1,41,sr) )     &
3035                                       * w(k,j,i) * rmask(j,i,sr)              &
3036                                       * rflags_invers(j,i,k+1)
3037                   ENDDO
3038                ENDDO
3039                sums_l(k,49,tn) = s1
3040             ENDDO
3041             !$acc end parallel loop
3042
3043          ENDIF
3044
3045          IF ( passive_scalar )  THEN
3046
3047             !$acc parallel loop gang present( hom, s, rflags_invers, rmask, sums_l, w ) create( s1 )
3048             DO  k = nzb, nzt_diff
3049                s1 = 0
3050                !$acc loop vector collapse( 2 ) reduction( +: s1 )
3051                DO  i = nxl, nxr
3052                   DO  j = nys, nyn
3053                      s1 = s1 + 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +      &
3054                                           s(k+1,j,i) - hom(k+1,1,117,sr) )    &
3055                                       * w(k,j,i) * rmask(j,i,sr)              &
3056                                       * rflags_invers(j,i,k+1)
3057                   ENDDO
3058                ENDDO
3059                sums_l(k,116,tn) = s1
3060             ENDDO
3061             !$acc end parallel loop
3062
3063          ENDIF
3064
3065       ENDIF
3066
3067
3068!
3069!--    Density at top follows Neumann condition
3070       IF ( ocean )  THEN
3071          !$acc parallel present( sums_l )
3072          sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
3073          sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
3074          !$acc end parallel
3075       ENDIF
3076
3077!
3078!--    Divergence of vertical flux of resolved scale energy and pressure
3079!--    fluctuations as well as flux of pressure fluctuation itself (68).
3080!--    First calculate the products, then the divergence.
3081!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
3082       IF ( hom(nzb+1,2,55,0) /= 0.0_wp  .OR.  hom(nzb+1,2,68,0) /= 0.0_wp )  THEN
3083
3084          STOP '+++ openACC porting for vertical flux div of resolved scale TKE in flow_statistics is still missing'
3085          sums_ll = 0.0_wp  ! local array
3086
3087          !$OMP DO
3088          DO  i = nxl, nxr
3089             DO  j = nys, nyn
3090                DO  k = nzb_s_inner(j,i)+1, nzt
3091
3092                   sums_ll(k,1) = sums_ll(k,1) + 0.5_wp * w(k,j,i) * (         &
3093                  ( 0.25_wp * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1) )  &
3094                            - 0.5_wp * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) ) )**2&
3095                + ( 0.25_wp * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i) )  &
3096                            - 0.5_wp * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) ) )**2&
3097                + w(k,j,i)**2                                        )
3098
3099                   sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i)             &
3100                                               * ( p(k,j,i) + p(k+1,j,i) )
3101
3102                ENDDO
3103             ENDDO
3104          ENDDO
3105          sums_ll(0,1)     = 0.0_wp    ! because w is zero at the bottom
3106          sums_ll(nzt+1,1) = 0.0_wp
3107          sums_ll(0,2)     = 0.0_wp
3108          sums_ll(nzt+1,2) = 0.0_wp
3109
3110          DO  k = nzb+1, nzt
3111             sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
3112             sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
3113             sums_l(k,68,tn) = sums_ll(k,2)
3114          ENDDO
3115          sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
3116          sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
3117          sums_l(nzb,68,tn) = 0.0_wp    ! because w* = 0 at nzb
3118
3119       ENDIF
3120
3121!
3122!--    Divergence of vertical flux of SGS TKE and the flux itself (69)
3123       IF ( hom(nzb+1,2,57,0) /= 0.0_wp  .OR.  hom(nzb+1,2,69,0) /= 0.0_wp )  THEN
3124
3125          STOP '+++ openACC porting for vertical flux div of SGS TKE in flow_statistics is still missing'
3126          !$OMP DO
3127          DO  i = nxl, nxr
3128             DO  j = nys, nyn
3129                DO  k = nzb_s_inner(j,i)+1, nzt
3130
3131                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5_wp * (              &
3132                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
3133                 - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
3134                                                                ) * ddzw(k)
3135
3136                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5_wp * (              &
3137                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
3138                                                                )
3139
3140                ENDDO
3141             ENDDO
3142          ENDDO
3143          sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
3144          sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
3145
3146       ENDIF
3147
3148!
3149!--    Horizontal heat fluxes (subgrid, resolved, total).
3150!--    Do it only, if profiles shall be plotted.
3151       IF ( hom(nzb+1,2,58,0) /= 0.0_wp ) THEN
3152
3153          STOP '+++ openACC porting for horizontal flux calculation in flow_statistics is still missing'
3154          !$OMP DO
3155          DO  i = nxl, nxr
3156             DO  j = nys, nyn
3157                DO  k = nzb_s_inner(j,i)+1, nzt
3158!
3159!--                Subgrid horizontal heat fluxes u"pt", v"pt"
3160                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5_wp *                &
3161                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
3162                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
3163                                                 * ddx * rmask(j,i,sr)
3164                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
3165                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
3166                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
3167                                                 * ddy * rmask(j,i,sr)
3168!
3169!--                Resolved horizontal heat fluxes u*pt*, v*pt*
3170                   sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
3171                                     ( u(k,j,i) - hom(k,1,1,sr) ) * 0.5_wp *   &
3172                                     ( pt(k,j,i-1) - hom(k,1,4,sr) +           &
3173                                       pt(k,j,i)   - hom(k,1,4,sr) )
3174                   pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
3175                                    pt(k,j,i)   - hom(k,1,4,sr) )
3176                   sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
3177                                     ( v(k,j,i) - hom(k,1,2,sr) ) * 0.5_wp *   &
3178                                     ( pt(k,j-1,i) - hom(k,1,4,sr) +           &
3179                                       pt(k,j,i)   - hom(k,1,4,sr) )
3180                ENDDO
3181             ENDDO
3182          ENDDO
3183!
3184!--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
3185          sums_l(nzb,58,tn) = 0.0_wp
3186          sums_l(nzb,59,tn) = 0.0_wp
3187          sums_l(nzb,60,tn) = 0.0_wp
3188          sums_l(nzb,61,tn) = 0.0_wp
3189          sums_l(nzb,62,tn) = 0.0_wp
3190          sums_l(nzb,63,tn) = 0.0_wp
3191
3192       ENDIF
3193
3194!
3195!--    Collect current large scale advection and subsidence tendencies for
3196!--    data output
3197       IF ( large_scale_forcing  .AND.  ( simulated_time > 0.0_wp ) )  THEN
3198!
3199!--       Interpolation in time of LSF_DATA
3200          nt = 1
3201          DO WHILE ( simulated_time - dt_3d > time_vert(nt) )
3202             nt = nt + 1
3203          ENDDO
3204          IF ( simulated_time - dt_3d /= time_vert(nt) )  THEN
3205            nt = nt - 1
3206          ENDIF
3207
3208          fac = ( simulated_time - dt_3d - time_vert(nt) )                     &
3209                / ( time_vert(nt+1)-time_vert(nt) )
3210
3211
3212          DO  k = nzb, nzt
3213             sums_ls_l(k,0) = td_lsa_lpt(k,nt)                                 &
3214                              + fac * ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )
3215             sums_ls_l(k,1) = td_lsa_q(k,nt)                                   &
3216                              + fac * ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) )
3217          ENDDO
3218
3219          sums_ls_l(nzt+1,0) = sums_ls_l(nzt,0)
3220          sums_ls_l(nzt+1,1) = sums_ls_l(nzt,1)
3221
3222          IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
3223
3224             DO  k = nzb, nzt
3225                sums_ls_l(k,2) = td_sub_lpt(k,nt) + fac *                      &
3226                                 ( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )
3227                sums_ls_l(k,3) = td_sub_q(k,nt) + fac *                        &
3228                                 ( td_sub_q(k,nt+1) - td_sub_q(k,nt) )
3229             ENDDO
3230
3231             sums_ls_l(nzt+1,2) = sums_ls_l(nzt,2)
3232             sums_ls_l(nzt+1,3) = sums_ls_l(nzt,3)
3233
3234          ENDIF
3235
3236       ENDIF
3237
3238
3239       IF ( land_surface )  THEN
3240          !$OMP DO
3241          DO  i = nxl, nxr
3242             DO  j =  nys, nyn
3243                DO  k = nzb_soil, nzt_soil
3244                   sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil(k,j,i)         &
3245                                      * rmask(j,i,sr)
3246                   sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil(k,j,i)         &
3247                                      * rmask(j,i,sr)
3248                ENDDO
3249             ENDDO
3250          ENDDO
3251       ENDIF
3252
3253
3254       IF ( radiation .AND. radiation_scheme == 'rrtmg' )  THEN
3255          !$OMP DO
3256          DO  i = nxl, nxr
3257             DO  j =  nys, nyn
3258                DO  k = nzb_s_inner(j,i)+1, nzt+1
3259                   sums_l(k,102,tn)  = sums_l(k,102,tn)  + rad_lw_in(k,j,i)    &
3260                                       * rmask(j,i,sr)
3261                   sums_l(k,103,tn)  = sums_l(k,103,tn)  + rad_lw_out(k,j,i)   &
3262                                       * rmask(j,i,sr)
3263                   sums_l(k,104,tn)  = sums_l(k,104,tn)  + rad_sw_in(k,j,i)    &
3264                                       * rmask(j,i,sr)
3265                   sums_l(k,105,tn)  = sums_l(k,105,tn)  + rad_sw_out(k,j,i)   &
3266                                       * rmask(j,i,sr)
3267#if defined ( __rrtmg )
3268                   sums_l(k,106,tn)  = sums_l(k,106,tn)  + rad_lw_cs_hr(k,j,i) &
3269                                       * rmask(j,i,sr)
3270                   sums_l(k,107,tn)  = sums_l(k,107,tn)  + rad_lw_hr(k,j,i)    &
3271                                       * rmask(j,i,sr)
3272                   sums_l(k,108,tn)  = sums_l(k,108,tn)  + rad_sw_cs_hr(k,j,i) &
3273                                       * rmask(j,i,sr)
3274                   sums_l(k,109,tn)  = sums_l(k,109,tn)  + rad_sw_hr(k,j,i)    &
3275                                       * rmask(j,i,sr)
3276#endif
3277                ENDDO
3278             ENDDO
3279          ENDDO
3280       ENDIF
3281
3282!
3283!--    Calculate the user-defined profiles
3284       CALL user_statistics( 'profiles', sr, tn )
3285       !$OMP END PARALLEL
3286
3287!
3288!--    Summation of thread sums
3289       IF ( threads_per_task > 1 )  THEN
3290          STOP '+++ openACC porting for threads_per_task > 1 in flow_statistics is still missing'
3291          DO  i = 1, threads_per_task-1
3292             sums_l(:,3,0)          = sums_l(:,3,0) + sums_l(:,3,i)
3293             sums_l(:,4:40,0)       = sums_l(:,4:40,0) + sums_l(:,4:40,i)
3294             sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
3295                                      sums_l(:,45:pr_palm,i)
3296             IF ( max_pr_user > 0 )  THEN
3297                sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
3298                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
3299                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
3300             ENDIF
3301          ENDDO
3302       ENDIF
3303
3304       !$acc update host( hom, sums, sums_l )
3305
3306#if defined( __parallel )
3307
3308!
3309!--    Compute total sum from local sums
3310       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
3311       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, &
3312                           MPI_SUM, comm2d, ierr )
3313       IF ( large_scale_forcing )  THEN
3314          CALL MPI_ALLREDUCE( sums_ls_l(nzb,2), sums(nzb,83), ngp_sums_ls,     &
3315                              MPI_REAL, MPI_SUM, comm2d, ierr )
3316       ENDIF
3317#else
3318       sums = sums_l(:,:,0)
3319       IF ( large_scale_forcing )  THEN
3320          sums(:,81:88) = sums_ls_l
3321       ENDIF
3322#endif
3323
3324!
3325!--    Final values are obtained by division by the total number of grid points
3326!--    used for summation. After that store profiles.
3327!--    Check, if statistical regions do contain at least one grid point at the
3328!--    respective k-level, otherwise division by zero will lead to undefined
3329!--    values, which may cause e.g. problems with NetCDF output
3330!--    Profiles:
3331       DO  k = nzb, nzt+1
3332          sums(k,3)             = sums(k,3)             / ngp_2dh(sr)
3333          sums(k,12:22)         = sums(k,12:22)         / ngp_2dh(sr)
3334          sums(k,30:32)         = sums(k,30:32)         / ngp_2dh(sr)
3335          sums(k,35:39)         = sums(k,35:39)         / ngp_2dh(sr)
3336          sums(k,45:53)         = sums(k,45:53)         / ngp_2dh(sr)
3337          sums(k,55:63)         = sums(k,55:63)         / ngp_2dh(sr)
3338          sums(k,81:88)         = sums(k,81:88)         / ngp_2dh(sr)
3339          sums(k,89:114)        = sums(k,89:114)        / ngp_2dh(sr)
3340          IF ( ngp_2dh_s_inner(k,sr) /= 0 )  THEN
3341             sums(k,8:11)          = sums(k,8:11)          / ngp_2dh_s_inner(k,sr)
3342             sums(k,23:29)         = sums(k,23:29)         / ngp_2dh_s_inner(k,sr)
3343             sums(k,33:34)         = sums(k,33:34)         / ngp_2dh_s_inner(k,sr)
3344             sums(k,40)            = sums(k,40)            / ngp_2dh_s_inner(k,sr)
3345             sums(k,54)            = sums(k,54)            / ngp_2dh_s_inner(k,sr)
3346             sums(k,64)            = sums(k,64)            / ngp_2dh_s_inner(k,sr)
3347             sums(k,70:80)         = sums(k,70:80)         / ngp_2dh_s_inner(k,sr)
3348             sums(k,115:pr_palm-2) = sums(k,115:pr_palm-2) / ngp_2dh_s_inner(k,sr)
3349          ENDIF
3350       ENDDO
3351
3352!--    u* and so on
3353!--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
3354!--    size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
3355!--    above the topography, they are being divided by ngp_2dh(sr)
3356       sums(nzb:nzb+3,pr_palm)    = sums(nzb:nzb+3,pr_palm)    / &
3357                                    ngp_2dh(sr)
3358       sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)       / &    ! qs
3359                                    ngp_2dh(sr)
3360!--    eges, e*
3361       sums(nzb+4:nzb+5,pr_palm)  = sums(nzb+4:nzb+5,pr_palm)  / &
3362                                    ngp_3d(sr)
3363!--    Old and new divergence
3364       sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
3365                                    ngp_3d_inner(sr)
3366
3367!--    User-defined profiles
3368       IF ( max_pr_user > 0 )  THEN
3369          DO  k = nzb, nzt+1
3370             IF ( ngp_2dh_s_inner(k,sr) /= 0 )  THEN
3371                sums(k,pr_palm+1:pr_palm+max_pr_user) = &
3372                                       sums(k,pr_palm+1:pr_palm+max_pr_user) / &
3373                                       ngp_2dh_s_inner(k,sr)
3374             ENDIF
3375          ENDDO
3376       ENDIF
3377
3378!
3379!--    Collect horizontal average in hom.
3380!--    Compute deduced averages (e.g. total heat flux)
3381       hom(:,1,3,sr)  = sums(:,3)      ! w
3382       hom(:,1,8,sr)  = sums(:,8)      ! e     profiles 5-7 are initial profiles
3383       hom(:,1,9,sr)  = sums(:,9)      ! km
3384       hom(:,1,10,sr) = sums(:,10)     ! kh
3385       hom(:,1,11,sr) = sums(:,11)     ! l
3386       hom(:,1,12,sr) = sums(:,12)     ! w"u"
3387       hom(:,1,13,sr) = sums(:,13)     ! w*u*
3388       hom(:,1,14,sr) = sums(:,14)     ! w"v"
3389       hom(:,1,15,sr) = sums(:,15)     ! w*v*
3390       hom(:,1,16,sr) = sums(:,16)     ! w"pt"
3391       hom(:,1,17,sr) = sums(:,17)     ! w*pt*
3392       hom(:,1,18,sr) = sums(:,16) + sums(:,17)    ! wpt
3393       hom(:,1,19,sr) = sums(:,12) + sums(:,13)    ! wu
3394       hom(:,1,20,sr) = sums(:,14) + sums(:,15)    ! wv
3395       hom(:,1,21,sr) = sums(:,21)     ! w*pt*BC
3396       hom(:,1,22,sr) = sums(:,16) + sums(:,21)    ! wptBC
3397                                       ! profile 24 is initial profile (sa)
3398                                       ! profiles 25-29 left empty for initial
3399                                       ! profiles
3400       hom(:,1,30,sr) = sums(:,30)     ! u*2
3401       hom(:,1,31,sr) = sums(:,31)     ! v*2
3402       hom(:,1,32,sr) = sums(:,32)     ! w*2
3403       hom(:,1,33,sr) = sums(:,33)     ! pt*2
3404       hom(:,1,34,sr) = sums(:,34)     ! e*
3405       hom(:,1,35,sr) = sums(:,35)     ! w*2pt*
3406       hom(:,1,36,sr) = sums(:,36)     ! w*pt*2
3407       hom(:,1,37,sr) = sums(:,37)     ! w*e*
3408       hom(:,1,38,sr) = sums(:,38)     ! w*3
3409       hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20_wp )**1.5_wp   ! Sw
3410       hom(:,1,40,sr) = sums(:,40)     ! p
3411       hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
3412       hom(:,1,46,sr) = sums(:,46)     ! w*vpt*
3413       hom(:,1,47,sr) = sums(:,45) + sums(:,46)    ! wvpt
3414       hom(:,1,48,sr) = sums(:,48)     ! w"q" (w"qv")
3415       hom(:,1,49,sr) = sums(:,49)     ! w*q* (w*qv*)
3416       hom(:,1,50,sr) = sums(:,48) + sums(:,49)    ! wq (wqv)
3417       hom(:,1,51,sr) = sums(:,51)     ! w"qv"
3418       hom(:,1,52,sr) = sums(:,52)     ! w*qv*
3419       hom(:,1,53,sr) = sums(:,52) + sums(:,51)    ! wq (wqv)
3420       hom(:,1,54,sr) = sums(:,54)     ! ql
3421       hom(:,1,55,sr) = sums(:,55)     ! w*u*u*/dz
3422       hom(:,1,56,sr) = sums(:,56)     ! w*p*/dz
3423       hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho )/dz
3424       hom(:,1,58,sr) = sums(:,58)     ! u"pt"
3425       hom(:,1,59,sr) = sums(:,59)     ! u*pt*
3426       hom(:,1,60,sr) = sums(:,58) + sums(:,59)    ! upt_t
3427       hom(:,1,61,sr) = sums(:,61)     ! v"pt"
3428       hom(:,1,62,sr) = sums(:,62)     ! v*pt*
3429       hom(:,1,63,sr) = sums(:,61) + sums(:,62)    ! vpt_t
3430       hom(:,1,64,sr) = sums(:,64)     ! rho
3431       hom(:,1,65,sr) = sums(:,65)     ! w"sa"
3432       hom(:,1,66,sr) = sums(:,66)     ! w*sa*
3433       hom(:,1,67,sr) = sums(:,65) + sums(:,66)    ! wsa
3434       hom(:,1,68,sr) = sums(:,68)     ! w*p*
3435       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho
3436       hom(:,1,70,sr) = sums(:,70)     ! q*2
3437       hom(:,1,71,sr) = sums(:,71)     ! prho
3438       hom(:,1,72,sr) = hyp * 1E-4_wp     ! hyp in dbar
3439       hom(:,1,73,sr) = sums(:,73)     ! nr
3440       hom(:,1,74,sr) = sums(:,74)     ! qr
3441       hom(:,1,75,sr) = sums(:,75)     ! qc
3442       hom(:,1,76,sr) = sums(:,76)     ! prr (precipitation rate)
3443                                       ! 77 is initial density profile
3444       hom(:,1,78,sr) = ug             ! ug
3445       hom(:,1,79,sr) = vg             ! vg
3446       hom(:,1,80,sr) = w_subs         ! w_subs
3447
3448       IF ( large_scale_forcing )  THEN
3449          hom(:,1,81,sr) = sums_ls_l(:,0)          ! td_lsa_lpt
3450          hom(:,1,82,sr) = sums_ls_l(:,1)          ! td_lsa_q
3451          IF ( use_subsidence_tendencies )  THEN
3452             hom(:,1,83,sr) = sums_ls_l(:,2)       ! td_sub_lpt
3453             hom(:,1,84,sr) = sums_ls_l(:,3)       ! td_sub_q
3454          ELSE
3455             hom(:,1,83,sr) = sums(:,83)           ! td_sub_lpt
3456             hom(:,1,84,sr) = sums(:,84)           ! td_sub_q
3457          ENDIF
3458          hom(:,1,85,sr) = sums(:,85)              ! td_nud_lpt
3459          hom(:,1,86,sr) = sums(:,86)              ! td_nud_q
3460          hom(:,1,87,sr) = sums(:,87)              ! td_nud_u
3461          hom(:,1,88,sr) = sums(:,88)              ! td_nud_v
3462       END IF
3463
3464       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
3465                                       ! u*, w'u', w'v', t* (in last profile)
3466
3467       IF ( max_pr_user > 0 )  THEN    ! user-defined profiles
3468          hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
3469                               sums(:,pr_palm+1:pr_palm+max_pr_user)
3470       ENDIF
3471
3472!
3473!--    Determine the boundary layer height using two different schemes.
3474!--    First scheme: Starting from the Earth's (Ocean's) surface, look for the
3475!--    first relative minimum (maximum) of the total heat flux.
3476!--    The corresponding height is assumed as the boundary layer height, if it
3477!--    is less than 1.5 times the height where the heat flux becomes negative
3478!--    (positive) for the first time.
3479       z_i(1) = 0.0_wp
3480       first = .TRUE.
3481
3482       IF ( ocean )  THEN
3483          DO  k = nzt, nzb+1, -1
3484             IF (  first  .AND.  hom(k,1,18,sr) < -1.0E-8_wp )  THEN
3485                first = .FALSE.
3486                height = zw(k)
3487             ENDIF
3488             IF ( hom(k,1,18,sr) < -1.0E-8_wp  .AND.                           &
3489                  hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
3490                IF ( zw(k) < 1.5_wp * height )  THEN
3491                   z_i(1) = zw(k)
3492                ELSE
3493                   z_i(1) = height
3494                ENDIF
3495                EXIT
3496             ENDIF
3497          ENDDO
3498       ELSE
3499          DO  k = nzb, nzt-1
3500             IF ( first  .AND.  hom(k,1,18,sr) < -1.0E-8_wp )  THEN
3501                first = .FALSE.
3502                height = zw(k)
3503             ENDIF
3504             IF ( hom(k,1,18,sr) < -1.0E-8_wp  .AND.                           &
3505                  hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
3506                IF ( zw(k) < 1.5_wp * height )  THEN
3507                   z_i(1) = zw(k)
3508                ELSE
3509                   z_i(1) = height
3510                ENDIF
3511                EXIT
3512             ENDIF
3513          ENDDO
3514       ENDIF
3515
3516!
3517!--    Second scheme: Gradient scheme from Sullivan et al. (1998), modified
3518!--    by Uhlenbrock(2006). The boundary layer height is the height with the
3519!--    maximal local temperature gradient: starting from the second (the last
3520!--    but one) vertical gridpoint, the local gradient must be at least
3521!--    0.2K/100m and greater than the next four gradients.
3522!--    WARNING: The threshold value of 0.2K/100m must be adjusted for the
3523!--             ocean case!
3524       z_i(2) = 0.0_wp
3525       DO  k = nzb+1, nzt+1
3526          dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
3527       ENDDO
3528       dptdz_threshold = 0.2_wp / 100.0_wp
3529
3530       IF ( ocean )  THEN
3531          DO  k = nzt+1, nzb+5, -1
3532             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
3533                  dptdz(k) > dptdz(k-1)  .AND.  dptdz(k) > dptdz(k-2)  .AND.  &
3534                  dptdz(k) > dptdz(k-3)  .AND.  dptdz(k) > dptdz(k-4) )  THEN
3535                z_i(2) = zw(k-1)
3536                EXIT
3537             ENDIF
3538          ENDDO
3539       ELSE
3540          DO  k = nzb+1, nzt-3
3541             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
3542                  dptdz(k) > dptdz(k+1)  .AND.  dptdz(k) > dptdz(k+2)  .AND.  &
3543                  dptdz(k) > dptdz(k+3)  .AND.  dptdz(k) > dptdz(k+4) )  THEN
3544                z_i(2) = zw(k-1)
3545                EXIT
3546             ENDIF
3547          ENDDO
3548       ENDIF
3549
3550       hom(nzb+6,1,pr_palm,sr) = z_i(1)
3551       hom(nzb+7,1,pr_palm,sr) = z_i(2)
3552
3553!
3554!--    Determine vertical index which is nearest to the mean surface level
3555!--    height of the respective statistic region
3556       DO  k = nzb, nzt
3557          IF ( zw(k) >= mean_surface_level_height(sr) )  THEN
3558             k_surface_level = k
3559             EXIT
3560          ENDIF
3561       ENDDO
3562
3563!
3564!--    Computation of both the characteristic vertical velocity and
3565!--    the characteristic convective boundary layer temperature.
3566!--    The inversion height entering into the equation is defined with respect
3567!--    to the mean surface level height of the respective statistic region.
3568!--    The horizontal average at surface level index + 1 is input for the
3569!--    average temperature.
3570       IF ( hom(nzb,1,18,sr) > 1.0E-8_wp  .AND.  z_i(1) /= 0.0_wp )  THEN
3571          hom(nzb+8,1,pr_palm,sr) = &
3572             ( g / hom(k_surface_level+1,1,4,sr) * hom(k_surface_level,1,18,sr)&
3573             * ABS( z_i(1) - mean_surface_level_height(sr) ) )**0.333333333_wp
3574       ELSE
3575          hom(nzb+8,1,pr_palm,sr)  = 0.0_wp
3576       ENDIF
3577
3578!
3579!--    Collect the time series quantities
3580       ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)     ! E
3581       ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)     ! E*
3582       ts_value(3,sr) = dt_3d
3583       ts_value(4,sr) = hom(nzb,1,pr_palm,sr)       ! u*
3584       ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr)     ! th*
3585       ts_value(6,sr) = u_max
3586       ts_value(7,sr) = v_max
3587       ts_value(8,sr) = w_max
3588       ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr)    ! new divergence
3589       ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr)    ! old Divergence
3590       ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr)    ! z_i(1)
3591       ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr)    ! z_i(2)
3592       ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr)    ! w*
3593       ts_value(14,sr) = hom(nzb,1,16,sr)           ! w'pt'   at k=0
3594       ts_value(15,sr) = hom(nzb+1,1,16,sr)         ! w'pt'   at k=1
3595       ts_value(16,sr) = hom(nzb+1,1,18,sr)         ! wpt     at k=1
3596       ts_value(17,sr) = hom(nzb,1,4,sr)            ! pt(0)
3597       ts_value(18,sr) = hom(nzb+1,1,4,sr)          ! pt(zp)
3598       ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr)    ! u'w'    at k=0
3599       ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)    ! v'w'    at k=0
3600       ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
3601
3602       IF ( .NOT. neutral )  THEN
3603          ts_value(22,sr) = hom(nzb,1,114,sr)          ! L
3604       ELSE
3605          ts_value(22,sr) = 1.0E10_wp
3606       ENDIF
3607
3608       ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
3609
3610!
3611!--    Collect land surface model timeseries
3612       IF ( land_surface )  THEN
3613          ts_value(dots_soil  ,sr) = hom(nzb,1,93,sr)           ! ghf_eb
3614          ts_value(dots_soil+1,sr) = hom(nzb,1,94,sr)           ! shf_eb
3615          ts_value(dots_soil+2,sr) = hom(nzb,1,95,sr)           ! qsws_eb
3616          ts_value(dots_soil+3,sr) = hom(nzb,1,96,sr)           ! qsws_liq_eb
3617          ts_value(dots_soil+4,sr) = hom(nzb,1,97,sr)           ! qsws_soil_eb
3618          ts_value(dots_soil+5,sr) = hom(nzb,1,98,sr)           ! qsws_veg_eb
3619          ts_value(dots_soil+6,sr) = hom(nzb,1,99,sr)           ! r_a
3620          ts_value(dots_soil+7,sr) = hom(nzb,1,100,sr)          ! r_s
3621       ENDIF
3622!
3623!--    Collect radiation model timeseries
3624       IF ( radiation )  THEN
3625          ts_value(dots_rad,sr)   = hom(nzb,1,101,sr)          ! rad_net
3626          ts_value(dots_rad+1,sr) = hom(nzb,1,102,sr)          ! rad_lw_in
3627          ts_value(dots_rad+2,sr) = hom(nzb,1,103,sr)          ! rad_lw_out
3628          ts_value(dots_rad+3,sr) = hom(nzb,1,104,sr)          ! rad_sw_in
3629          ts_value(dots_rad+4,sr) = hom(nzb,1,105,sr)          ! rad_sw_out
3630
3631          IF ( radiation_scheme == 'rrtmg' )  THEN
3632             ts_value(dots_rad+5,sr) = hom(nzb,1,106,sr)          ! rrtm_aldif
3633             ts_value(dots_rad+6,sr) = hom(nzb,1,107,sr)          ! rrtm_aldir
3634             ts_value(dots_rad+7,sr) = hom(nzb,1,108,sr)          ! rrtm_asdif
3635             ts_value(dots_rad+8,sr) = hom(nzb,1,109,sr)          ! rrtm_asdir
3636          ENDIF
3637
3638       ENDIF
3639
3640!
3641!--    Calculate additional statistics provided by the user interface
3642       CALL user_statistics( 'time_series', sr, 0 )
3643
3644    ENDDO    ! loop of the subregions
3645
3646    !$acc end data
3647
3648!
3649!-- If required, sum up horizontal averages for subsequent time averaging
3650!-- Do not sum, if flow statistics is called before the first initial time step.
3651    IF ( do_sum  .AND.  simulated_time /= 0.0_wp )  THEN
3652       IF ( average_count_pr == 0 )  hom_sum = 0.0_wp
3653       hom_sum = hom_sum + hom(:,1,:,:)
3654       average_count_pr = average_count_pr + 1
3655       do_sum = .FALSE.
3656    ENDIF
3657
3658!
3659!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
3660!-- This flag is reset after each time step in time_integration.
3661    flow_statistics_called = .TRUE.
3662
3663    CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
3664
3665
3666 END SUBROUTINE flow_statistics
3667#endif
Note: See TracBrowser for help on using the repository browser.