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

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

bugfixes for calculating run control quantities, bugfix for calculating pressure with fft-method in case of Neumann conditions both at bottom and top, steering of pres modified, ocean mode now uses initial density profile as reference in the buoyancy term

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