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

Last change on this file since 1655 was 1655, checked in by raasch, 6 years ago

last commit documented

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