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

Last change on this file since 1710 was 1710, checked in by maronga, 6 years ago

last commit documented

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