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

Last change on this file since 1783 was 1783, checked in by raasch, 8 years ago

NetCDF routines modularized; new parameter netcdf_deflate; further changes in the pmc

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