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

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

FORTRAN error removed from openacc branch

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