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

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

last commit documented

  • Property svn:keywords set to Id
File size: 142.4 KB
Line 
1#if ! defined( __openacc )
2 SUBROUTINE flow_statistics
3
4!--------------------------------------------------------------------------------!
5! This file is part of PALM.
6!
7! PALM is free software: you can redistribute it and/or modify it under the terms
8! of the GNU General Public License as published by the Free Software Foundation,
9! either version 3 of the License, or (at your option) any later version.
10!
11! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
12! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
13! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
14!
15! You should have received a copy of the GNU General Public License along with
16! PALM. If not, see <http://www.gnu.org/licenses/>.
17!
18! Copyright 1997-2014 Leibniz Universitaet Hannover
19!--------------------------------------------------------------------------------!
20!
21! Current revisions:
22! -----------------
23!
24!
25! Former revisions:
26! -----------------
27! $Id: flow_statistics.f90 1586 2015-04-30 07:11:33Z maronga $
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 pegrid
1615   
1616    USE statistics
1617
1618    IMPLICIT NONE
1619
1620    INTEGER(iwp) ::  i                   !:
1621    INTEGER(iwp) ::  j                   !:
1622    INTEGER(iwp) ::  k                   !:
1623    INTEGER(iwp) ::  nt                  !:
1624    INTEGER(iwp) ::  omp_get_thread_num  !:
1625    INTEGER(iwp) ::  sr                  !:
1626    INTEGER(iwp) ::  tn                  !:
1627   
1628    LOGICAL ::  first  !:
1629   
1630    REAL(wp) ::  dptdz_threshold  !:
1631    REAL(wp) ::  fac              !:
1632    REAL(wp) ::  height           !:
1633    REAL(wp) ::  pts              !:
1634    REAL(wp) ::  sums_l_eper      !:
1635    REAL(wp) ::  sums_l_etot      !:
1636    REAL(wp) ::  s1               !:
1637    REAL(wp) ::  s2               !:
1638    REAL(wp) ::  s3               !:
1639    REAL(wp) ::  s4               !:
1640    REAL(wp) ::  s5               !:
1641    REAL(wp) ::  s6               !:
1642    REAL(wp) ::  s7               !:
1643    REAL(wp) ::  ust              !:
1644    REAL(wp) ::  ust2             !:
1645    REAL(wp) ::  u2               !:
1646    REAL(wp) ::  vst              !:
1647    REAL(wp) ::  vst2             !:
1648    REAL(wp) ::  v2               !:
1649    REAL(wp) ::  w2               !:
1650    REAL(wp) ::  z_i(2)           !:
1651
1652    REAL(wp) ::  dptdz(nzb+1:nzt+1)    !:
1653    REAL(wp) ::  sums_ll(nzb:nzt+1,2)  !:
1654
1655    CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
1656
1657!
1658!-- To be on the safe side, check whether flow_statistics has already been
1659!-- called once after the current time step
1660    IF ( flow_statistics_called )  THEN
1661
1662       message_string = 'flow_statistics is called two times within one ' // &
1663                        'timestep'
1664       CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
1665
1666    ENDIF
1667
1668    !$acc data create( sums, sums_l )
1669    !$acc update device( hom )
1670
1671!
1672!-- Compute statistics for each (sub-)region
1673    DO  sr = 0, statistic_regions
1674
1675!
1676!--    Initialize (local) summation array
1677       sums_l = 0.0_wp
1678
1679!
1680!--    Store sums that have been computed in other subroutines in summation
1681!--    array
1682       sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
1683!--    WARNING: next line still has to be adjusted for OpenMP
1684       sums_l(:,21,0) = sums_wsts_bc_l(:,sr)  ! heat flux from advec_s_bc
1685       sums_l(nzb+9,pr_palm,0)  = sums_divold_l(sr)  ! old divergence from pres
1686       sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
1687
1688!
1689!--    When calcuating horizontally-averaged total (resolved- plus subgrid-
1690!--    scale) vertical fluxes and velocity variances by using commonly-
1691!--    applied Reynolds-based methods ( e.g. <w'pt'> = (w-<w>)*(pt-<pt>) )
1692!--    in combination with the 5th order advection scheme, pronounced
1693!--    artificial kinks could be observed in the vertical profiles near the
1694!--    surface. Please note: these kinks were not related to the model truth,
1695!--    i.e. these kinks are just related to an evaluation problem.   
1696!--    In order avoid these kinks, vertical fluxes and horizontal as well
1697!--    vertical velocity variances are calculated directly within the advection
1698!--    routines, according to the numerical discretization, to evaluate the
1699!--    statistical quantities as they will appear within the prognostic
1700!--    equations.
1701!--    Copy the turbulent quantities, evaluated in the advection routines to
1702!--    the local array sums_l() for further computations.
1703       IF ( ws_scheme_mom .AND. sr == 0 )  THEN
1704
1705!
1706!--       According to the Neumann bc for the horizontal velocity components,
1707!--       the corresponding fluxes has to satisfiy the same bc.
1708          IF ( ocean )  THEN
1709             sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
1710             sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)
1711          ENDIF
1712
1713          DO  i = 0, threads_per_task-1
1714!
1715!--          Swap the turbulent quantities evaluated in advec_ws.
1716             sums_l(:,13,i) = sums_wsus_ws_l(:,i)       ! w*u*
1717             sums_l(:,15,i) = sums_wsvs_ws_l(:,i)       ! w*v*
1718             sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2
1719             sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2
1720             sums_l(:,32,i) = sums_ws2_ws_l(:,i)        ! w*2
1721             sums_l(:,34,i) = sums_l(:,34,i) + 0.5_wp *                        &
1722                              ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +      &
1723                                sums_ws2_ws_l(:,i) )    ! e*
1724             DO  k = nzb, nzt
1725                sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5_wp * ( &
1726                                                      sums_us2_ws_l(k,i) +     &
1727                                                      sums_vs2_ws_l(k,i) +     &
1728                                                      sums_ws2_ws_l(k,i)     )
1729             ENDDO
1730          ENDDO
1731
1732       ENDIF
1733
1734       IF ( ws_scheme_sca .AND. sr == 0 )  THEN
1735
1736          DO  i = 0, threads_per_task-1
1737             sums_l(:,17,i) = sums_wspts_ws_l(:,i)      ! w*pt* from advec_s_ws
1738             IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
1739             IF ( humidity .OR. passive_scalar ) sums_l(:,49,i) =              &
1740                                                   sums_wsqs_ws_l(:,i) !w*q*
1741          ENDDO
1742
1743       ENDIF
1744!
1745!--    Horizontally averaged profiles of horizontal velocities and temperature.
1746!--    They must have been computed before, because they are already required
1747!--    for other horizontal averages.
1748       tn = 0
1749
1750       !$OMP PARALLEL PRIVATE( i, j, k, tn )
1751#if defined( __intel_openmp_bug )
1752       tn = omp_get_thread_num()
1753#else
1754!$     tn = omp_get_thread_num()
1755#endif
1756
1757       !$acc update device( sums_l )
1758
1759       !$OMP DO
1760       !$acc parallel loop gang present( pt, rflags_invers, rmask, sums_l, u, v ) create( s1, s2, s3 )
1761       DO  k = nzb, nzt+1
1762          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
1763          DO  i = nxl, nxr
1764             DO  j =  nys, nyn
1765!
1766!--             k+1 is used in rflags since rflags is set 0 at surface points
1767                s1 = s1 + u(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1768                s2 = s2 + v(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1769                s3 = s3 + pt(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1770             ENDDO
1771          ENDDO
1772          sums_l(k,1,tn) = s1
1773          sums_l(k,2,tn) = s2
1774          sums_l(k,4,tn) = s3
1775       ENDDO
1776       !$acc end parallel loop
1777
1778!
1779!--    Horizontally averaged profile of salinity
1780       IF ( ocean )  THEN
1781          !$OMP DO
1782          !$acc parallel loop gang present( rflags_invers, rmask, sums_l, sa ) create( s1 )
1783          DO  k = nzb, nzt+1
1784             !$acc loop vector collapse( 2 ) reduction( +: s1 )
1785             DO  i = nxl, nxr
1786                DO  j =  nys, nyn
1787                   s1 = s1 + sa(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1788                ENDDO
1789             ENDDO
1790             sums_l(k,23,tn) = s1
1791          ENDDO
1792          !$acc end parallel loop
1793       ENDIF
1794
1795!
1796!--    Horizontally averaged profiles of virtual potential temperature,
1797!--    total water content, specific humidity and liquid water potential
1798!--    temperature
1799       IF ( humidity )  THEN
1800
1801          !$OMP DO
1802          !$acc parallel loop gang present( q, rflags_invers, rmask, sums_l, vpt ) create( s1, s2 )
1803          DO  k = nzb, nzt+1
1804             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
1805             DO  i = nxl, nxr
1806                DO  j =  nys, nyn
1807                   s1 = s1 + q(k,j,i)   * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1808                   s2 = s2 + vpt(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1809                ENDDO
1810             ENDDO
1811             sums_l(k,41,tn) = s1
1812             sums_l(k,44,tn) = s2
1813          ENDDO
1814          !$acc end parallel loop
1815
1816          IF ( cloud_physics )  THEN
1817             !$OMP DO
1818             !$acc parallel loop gang present( pt, q, ql, rflags_invers, rmask, sums_l ) create( s1, s2 )
1819             DO  k = nzb, nzt+1
1820                !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
1821                DO  i = nxl, nxr
1822                   DO  j =  nys, nyn
1823                      s1 = s1 + ( q(k,j,i) - ql(k,j,i) ) * &
1824                                rmask(j,i,sr) * rflags_invers(j,i,k+1)
1825                      s2 = s2 + ( pt(k,j,i) + l_d_cp*pt_d_t(k) * ql(k,j,i) ) * &
1826                                rmask(j,i,sr) * rflags_invers(j,i,k+1)
1827                   ENDDO
1828                ENDDO
1829                sums_l(k,42,tn) = s1
1830                sums_l(k,43,tn) = s2
1831             ENDDO
1832             !$acc end parallel loop
1833          ENDIF
1834       ENDIF
1835
1836!
1837!--    Horizontally averaged profiles of passive scalar
1838       IF ( passive_scalar )  THEN
1839          !$OMP DO
1840          !$acc parallel loop gang present( q, rflags_invers, rmask, sums_l ) create( s1 )
1841          DO  k = nzb, nzt+1
1842             !$acc loop vector collapse( 2 ) reduction( +: s1 )
1843             DO  i = nxl, nxr
1844                DO  j =  nys, nyn
1845                   s1 = s1 + q(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1846                ENDDO
1847             ENDDO
1848             sums_l(k,41,tn) = s1
1849          ENDDO
1850          !$acc end parallel loop
1851       ENDIF
1852       !$OMP END PARALLEL
1853
1854!
1855!--    Summation of thread sums
1856       IF ( threads_per_task > 1 )  THEN
1857          DO  i = 1, threads_per_task-1
1858             !$acc parallel present( sums_l )
1859             sums_l(:,1,0) = sums_l(:,1,0) + sums_l(:,1,i)
1860             sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i)
1861             sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i)
1862             !$acc end parallel
1863             IF ( ocean )  THEN
1864                !$acc parallel present( sums_l )
1865                sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)
1866                !$acc end parallel
1867             ENDIF
1868             IF ( humidity )  THEN
1869                !$acc parallel present( sums_l )
1870                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
1871                sums_l(:,44,0) = sums_l(:,44,0) + sums_l(:,44,i)
1872                !$acc end parallel
1873                IF ( cloud_physics )  THEN
1874                   !$acc parallel present( sums_l )
1875                   sums_l(:,42,0) = sums_l(:,42,0) + sums_l(:,42,i)
1876                   sums_l(:,43,0) = sums_l(:,43,0) + sums_l(:,43,i)
1877                   !$acc end parallel
1878                ENDIF
1879             ENDIF
1880             IF ( passive_scalar )  THEN
1881                !$acc parallel present( sums_l )
1882                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
1883                !$acc end parallel
1884             ENDIF
1885          ENDDO
1886       ENDIF
1887
1888#if defined( __parallel )
1889!
1890!--    Compute total sum from local sums
1891       !$acc update host( sums_l )
1892       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1893       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL,  &
1894                           MPI_SUM, comm2d, ierr )
1895       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1896       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL,  &
1897                           MPI_SUM, comm2d, ierr )
1898       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1899       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL,  &
1900                           MPI_SUM, comm2d, ierr )
1901       IF ( ocean )  THEN
1902          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1903          CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb,       &
1904                              MPI_REAL, MPI_SUM, comm2d, ierr )
1905       ENDIF
1906       IF ( humidity ) THEN
1907          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1908          CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb,       &
1909                              MPI_REAL, MPI_SUM, comm2d, ierr )
1910          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1911          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb,       &
1912                              MPI_REAL, MPI_SUM, comm2d, ierr )
1913          IF ( cloud_physics ) THEN
1914             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1915             CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb,    &
1916                                 MPI_REAL, MPI_SUM, comm2d, ierr )
1917             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1918             CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb,    &
1919                                 MPI_REAL, MPI_SUM, comm2d, ierr )
1920          ENDIF
1921       ENDIF
1922
1923       IF ( passive_scalar )  THEN
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       ENDIF
1928       !$acc update device( sums )
1929#else
1930       !$acc parallel present( sums, sums_l )
1931       sums(:,1) = sums_l(:,1,0)
1932       sums(:,2) = sums_l(:,2,0)
1933       sums(:,4) = sums_l(:,4,0)
1934       !$acc end parallel
1935       IF ( ocean )  THEN
1936          !$acc parallel present( sums, sums_l )
1937          sums(:,23) = sums_l(:,23,0)
1938          !$acc end parallel
1939       ENDIF
1940       IF ( humidity )  THEN
1941          !$acc parallel present( sums, sums_l )
1942          sums(:,44) = sums_l(:,44,0)
1943          sums(:,41) = sums_l(:,41,0)
1944          !$acc end parallel
1945          IF ( cloud_physics )  THEN
1946             !$acc parallel present( sums, sums_l )
1947             sums(:,42) = sums_l(:,42,0)
1948             sums(:,43) = sums_l(:,43,0)
1949             !$acc end parallel
1950          ENDIF
1951       ENDIF
1952       IF ( passive_scalar )  THEN
1953          !$acc parallel present( sums, sums_l )
1954          sums(:,41) = sums_l(:,41,0)
1955          !$acc end parallel
1956       ENDIF
1957#endif
1958
1959!
1960!--    Final values are obtained by division by the total number of grid points
1961!--    used for summation. After that store profiles.
1962       !$acc parallel present( hom, ngp_2dh, ngp_2dh_s_inner, sums )
1963       sums(:,1) = sums(:,1) / ngp_2dh(sr)
1964       sums(:,2) = sums(:,2) / ngp_2dh(sr)
1965       sums(:,4) = sums(:,4) / ngp_2dh_s_inner(:,sr)
1966       hom(:,1,1,sr) = sums(:,1)             ! u
1967       hom(:,1,2,sr) = sums(:,2)             ! v
1968       hom(:,1,4,sr) = sums(:,4)             ! pt
1969       !$acc end parallel
1970
1971!
1972!--    Salinity
1973       IF ( ocean )  THEN
1974          !$acc parallel present( hom, ngp_2dh_s_inner, sums )
1975          sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr)
1976          hom(:,1,23,sr) = sums(:,23)             ! sa
1977          !$acc end parallel
1978       ENDIF
1979
1980!
1981!--    Humidity and cloud parameters
1982       IF ( humidity ) THEN
1983          !$acc parallel present( hom, ngp_2dh_s_inner, sums )
1984          sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr)
1985          sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
1986          hom(:,1,44,sr) = sums(:,44)                ! vpt
1987          hom(:,1,41,sr) = sums(:,41)                ! qv (q)
1988          !$acc end parallel
1989          IF ( cloud_physics ) THEN
1990             !$acc parallel present( hom, ngp_2dh_s_inner, sums )
1991             sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr)
1992             sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr)
1993             hom(:,1,42,sr) = sums(:,42)             ! qv
1994             hom(:,1,43,sr) = sums(:,43)             ! pt
1995             !$acc end parallel
1996          ENDIF
1997       ENDIF
1998
1999!
2000!--    Passive scalar
2001       IF ( passive_scalar )  THEN
2002          !$acc parallel present( hom, ngp_2dh_s_inner, sums )
2003          sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
2004          hom(:,1,41,sr) = sums(:,41)                ! s (q)
2005          !$acc end parallel
2006       ENDIF
2007
2008!
2009!--    Horizontally averaged profiles of the remaining prognostic variables,
2010!--    variances, the total and the perturbation energy (single values in last
2011!--    column of sums_l) and some diagnostic quantities.
2012!--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
2013!--    ----  speaking the following k-loop would have to be split up and
2014!--          rearranged according to the staggered grid.
2015!--          However, this implies no error since staggered velocity components
2016!--          are zero at the walls and inside buildings.
2017       tn = 0
2018#if defined( __intel_openmp_bug )
2019       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, &
2020       !$OMP                    tn, ust, ust2, u2, vst, vst2, v2, w2 )
2021       tn = omp_get_thread_num()
2022#else
2023       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2, w2 )
2024!$     tn = omp_get_thread_num()
2025#endif
2026       !$OMP DO
2027       !$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 )
2028       DO  k = nzb, nzt+1
2029          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5, s6, s7 )
2030          DO  i = nxl, nxr
2031             DO  j =  nys, nyn
2032!
2033!--             Prognostic and diagnostic variables
2034                s1 = s1 + w(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2035                s2 = s2 + e(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2036                s3 = s3 + km(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2037                s4 = s4 + kh(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2038                s5 = s5 + p(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2039                s6 = s6 + ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr) * &
2040                          rflags_invers(j,i,k+1)
2041!
2042!--             Higher moments
2043!--             (Computation of the skewness of w further below)
2044                s7 = s7 + w(k,j,i)**3 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2045             ENDDO
2046          ENDDO
2047          sums_l(k,3,tn)  = s1
2048          sums_l(k,8,tn)  = s2
2049          sums_l(k,9,tn)  = s3
2050          sums_l(k,10,tn) = s4
2051          sums_l(k,40,tn) = s5
2052          sums_l(k,33,tn) = s6
2053          sums_l(k,38,tn) = s7
2054       ENDDO
2055       !$acc end parallel loop
2056
2057       IF ( humidity )  THEN
2058          !$OMP DO
2059          !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l ) create( s1 )
2060          DO  k = nzb, nzt+1
2061             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2062             DO  i = nxl, nxr
2063                DO  j =  nys, nyn
2064                   s1 = s1 + ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr) * &
2065                             rflags_invers(j,i,k+1)
2066                ENDDO
2067             ENDDO
2068             sums_l(k,70,tn) = s1
2069          ENDDO
2070          !$acc end parallel loop
2071       ENDIF
2072
2073!
2074!--    Total and perturbation energy for the total domain (being
2075!--    collected in the last column of sums_l).
2076       !$OMP DO
2077       !$acc parallel loop collapse(3) present( rflags_invers, rmask, u, v, w ) reduction(+:s1)
2078       DO  i = nxl, nxr
2079          DO  j =  nys, nyn
2080             DO  k = nzb, nzt+1
2081                s1 = s1 + 0.5_wp *                                             & 
2082                          ( u(k,j,i)**2 + v(k,j,i)**2 + w(k,j,i)**2 ) *        &
2083                          rmask(j,i,sr) * rflags_invers(j,i,k+1)
2084             ENDDO
2085          ENDDO
2086       ENDDO
2087       !$acc end parallel loop
2088       !$acc parallel present( sums_l )
2089       sums_l(nzb+4,pr_palm,tn) = s1
2090       !$acc end parallel
2091
2092       !$OMP DO
2093       !$acc parallel present( rmask, sums_l, us, usws, vsws, ts ) create( s1, s2, s3, s4 )
2094       !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4 )
2095       DO  i = nxl, nxr
2096          DO  j =  nys, nyn
2097!
2098!--          2D-arrays (being collected in the last column of sums_l)
2099             s1 = s1 + us(j,i)   * rmask(j,i,sr)
2100             s2 = s2 + usws(j,i) * rmask(j,i,sr)
2101             s3 = s3 + vsws(j,i) * rmask(j,i,sr)
2102             s4 = s4 + ts(j,i)   * rmask(j,i,sr)
2103          ENDDO
2104       ENDDO
2105       sums_l(nzb,pr_palm,tn)   = s1
2106       sums_l(nzb+1,pr_palm,tn) = s2
2107       sums_l(nzb+2,pr_palm,tn) = s3
2108       sums_l(nzb+3,pr_palm,tn) = s4
2109       !$acc end parallel
2110
2111       IF ( humidity )  THEN
2112          !$acc parallel present( qs, rmask, sums_l ) create( s1 )
2113          !$acc loop vector collapse( 2 ) reduction( +: s1 )
2114          DO  i = nxl, nxr
2115             DO  j =  nys, nyn
2116                s1 = s1 + qs(j,i) * rmask(j,i,sr)
2117             ENDDO
2118          ENDDO
2119          sums_l(nzb+12,pr_palm,tn) = s1
2120          !$acc end parallel
2121       ENDIF
2122
2123!
2124!--    Computation of statistics when ws-scheme is not used. Else these
2125!--    quantities are evaluated in the advection routines.
2126       IF ( .NOT. ws_scheme_mom  .OR.  sr /= 0 )  THEN
2127
2128          !$OMP DO
2129          !$acc parallel loop gang present( u, v, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3, s4, ust2, vst2, w2 )
2130          DO  k = nzb, nzt+1
2131             !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4 )
2132             DO  i = nxl, nxr
2133                DO  j =  nys, nyn
2134                   ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
2135                   vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
2136                   w2   = w(k,j,i)**2
2137
2138                   s1 = s1 + ust2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2139                   s2 = s2 + vst2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2140                   s3 = s3 + w2   * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2141!
2142!--                Perturbation energy
2143                   s4 = s4 + 0.5_wp * ( ust2 + vst2 + w2 ) * rmask(j,i,sr) *   &
2144                             rflags_invers(j,i,k+1)
2145                ENDDO
2146             ENDDO
2147             sums_l(k,30,tn) = s1
2148             sums_l(k,31,tn) = s2
2149             sums_l(k,32,tn) = s3
2150             sums_l(k,34,tn) = s4
2151          ENDDO
2152          !$acc end parallel loop
2153!
2154!--       Total perturbation TKE
2155          !$OMP DO
2156          !$acc parallel present( sums_l ) create( s1 )
2157          !$acc loop reduction( +: s1 )
2158          DO  k = nzb, nzt+1
2159             s1 = s1 + sums_l(k,34,tn)
2160          ENDDO
2161          sums_l(nzb+5,pr_palm,tn) = s1
2162          !$acc end parallel
2163
2164       ENDIF
2165
2166!
2167!--    Horizontally averaged profiles of the vertical fluxes
2168
2169!
2170!--    Subgridscale fluxes.
2171!--    WARNING: If a Prandtl-layer is used (k=nzb for flat terrain), the fluxes
2172!--    -------  should be calculated there in a different way. This is done
2173!--             in the next loop further below, where results from this loop are
2174!--             overwritten. However, THIS WORKS IN CASE OF FLAT TERRAIN ONLY!
2175!--             The non-flat case still has to be handled.
2176!--    NOTE: for simplicity, nzb_s_inner is used below, although
2177!--    ----  strictly speaking the following k-loop would have to be
2178!--          split up according to the staggered grid.
2179!--          However, this implies no error since staggered velocity
2180!--          components are zero at the walls and inside buildings.
2181       !$OMP DO
2182       !$acc parallel loop gang present( ddzu, kh, km, pt, u, v, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3 )
2183       DO  k = nzb, nzt_diff
2184          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
2185          DO  i = nxl, nxr
2186             DO  j = nys, nyn
2187
2188!
2189!--             Momentum flux w"u"
2190                s1 = s1 - 0.25_wp * (                                          &
2191                               km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
2192                                                           ) * (               &
2193                                   ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
2194                                 + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
2195                                                               )               &
2196                               * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2197!
2198!--             Momentum flux w"v"
2199                s2 = s2 - 0.25_wp * (                                          &
2200                               km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
2201                                                           ) * (               &
2202                                   ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
2203                                 + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
2204                                                               )               &
2205                               * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2206!
2207!--             Heat flux w"pt"
2208                s3 = s3 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )                 &
2209                                 * ( pt(k+1,j,i) - pt(k,j,i) )                 &
2210                                 * ddzu(k+1) * rmask(j,i,sr)                   &
2211                                 * rflags_invers(j,i,k+1)
2212             ENDDO
2213          ENDDO
2214          sums_l(k,12,tn) = s1
2215          sums_l(k,14,tn) = s2
2216          sums_l(k,16,tn) = s3
2217       ENDDO
2218       !$acc end parallel loop
2219
2220!
2221!--    Salinity flux w"sa"
2222       IF ( ocean )  THEN
2223          !$acc parallel loop gang present( ddzu, kh, sa, rflags_invers, rmask, sums_l ) create( s1 )
2224          DO  k = nzb, nzt_diff
2225             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2226             DO  i = nxl, nxr
2227                DO  j = nys, nyn
2228                   s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
2229                                    * ( sa(k+1,j,i) - sa(k,j,i) )              &
2230                                    * ddzu(k+1) * rmask(j,i,sr)                & 
2231                                    * rflags_invers(j,i,k+1)
2232                ENDDO
2233             ENDDO
2234             sums_l(k,65,tn) = s1
2235          ENDDO
2236          !$acc end parallel loop
2237       ENDIF
2238
2239!
2240!--    Buoyancy flux, water flux (humidity flux) w"q"
2241       IF ( humidity ) THEN
2242
2243          !$acc parallel loop gang present( ddzu, kh, q, vpt, rflags_invers, rmask, sums_l ) create( s1, s2 )
2244          DO  k = nzb, nzt_diff
2245             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2246             DO  i = nxl, nxr
2247                DO  j = nys, nyn
2248                   s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
2249                                    * ( vpt(k+1,j,i) - vpt(k,j,i) )            &
2250                                    * ddzu(k+1) * rmask(j,i,sr)                &
2251                                    * rflags_invers(j,i,k+1)
2252                   s2 = s2 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
2253                                    * ( q(k+1,j,i) - q(k,j,i) )                &
2254                                    * ddzu(k+1) * rmask(j,i,sr)                &
2255                                    * rflags_invers(j,i,k+1)
2256                ENDDO
2257             ENDDO
2258             sums_l(k,45,tn) = s1
2259             sums_l(k,48,tn) = s2
2260          ENDDO
2261          !$acc end parallel loop
2262
2263          IF ( cloud_physics ) THEN
2264
2265             !$acc parallel loop gang present( ddzu, kh, q, ql, rflags_invers, rmask, sums_l ) create( s1 )
2266             DO  k = nzb, nzt_diff
2267                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2268                DO  i = nxl, nxr
2269                   DO  j = nys, nyn
2270                      s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )           &
2271                                       *  ( ( q(k+1,j,i) - ql(k+1,j,i) )       &
2272                                          - ( q(k,j,i) - ql(k,j,i) ) )         &
2273                                       * ddzu(k+1) * rmask(j,i,sr)             & 
2274                                       * rflags_invers(j,i,k+1)
2275                   ENDDO
2276                ENDDO
2277                sums_l(k,51,tn) = s1
2278             ENDDO
2279             !$acc end parallel loop
2280
2281          ENDIF
2282
2283       ENDIF
2284!
2285!--    Passive scalar flux
2286       IF ( passive_scalar )  THEN
2287
2288          !$acc parallel loop gang present( ddzu, kh, q, rflags_invers, rmask, sums_l ) create( s1 )
2289          DO  k = nzb, nzt_diff
2290             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2291             DO  i = nxl, nxr
2292                DO  j = nys, nyn
2293                   s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
2294                                    * ( q(k+1,j,i) - q(k,j,i) )                &
2295                                    * ddzu(k+1) * rmask(j,i,sr)                &
2296                                    * rflags_invers(j,i,k+1)
2297                ENDDO
2298             ENDDO
2299             sums_l(k,48,tn) = s1
2300          ENDDO
2301          !$acc end parallel loop
2302
2303       ENDIF
2304
2305       IF ( use_surface_fluxes )  THEN
2306
2307          !$OMP DO
2308          !$acc parallel present( rmask, shf, sums_l, usws, vsws ) create( s1, s2, s3, s4, s5 )
2309          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5 )
2310          DO  i = nxl, nxr
2311             DO  j =  nys, nyn
2312!
2313!--             Subgridscale fluxes in the Prandtl layer
2314                s1 = s1 + usws(j,i) * rmask(j,i,sr)     ! w"u"
2315                s2 = s2 + vsws(j,i) * rmask(j,i,sr)     ! w"v"
2316                s3 = s3 + shf(j,i)  * rmask(j,i,sr)     ! w"pt"
2317                s4 = s4 + 0.0_wp * rmask(j,i,sr)        ! u"pt"
2318                s5 = s5 + 0.0_wp * rmask(j,i,sr)        ! v"pt"
2319             ENDDO
2320          ENDDO
2321          sums_l(nzb,12,tn) = s1
2322          sums_l(nzb,14,tn) = s2
2323          sums_l(nzb,16,tn) = s3
2324          sums_l(nzb,58,tn) = s4
2325          sums_l(nzb,61,tn) = s5
2326          !$acc end parallel
2327
2328          IF ( ocean )  THEN
2329
2330             !$OMP DO
2331             !$acc parallel present( rmask, saswsb, sums_l ) create( s1 )
2332             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2333             DO  i = nxl, nxr
2334                DO  j =  nys, nyn
2335                   s1 = s1 + saswsb(j,i) * rmask(j,i,sr)  ! w"sa"
2336                ENDDO
2337             ENDDO
2338             sums_l(nzb,65,tn) = s1
2339             !$acc end parallel
2340
2341          ENDIF
2342
2343          IF ( humidity )  THEN
2344
2345             !$OMP DO
2346             !$acc parallel present( pt, q, qsws, rmask, shf, sums_l ) create( s1, s2 )
2347             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2348             DO  i = nxl, nxr
2349                DO  j =  nys, nyn
2350                   s1 = s1 + qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
2351                   s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) * shf(j,i)    &
2352                               + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )
2353                ENDDO
2354             ENDDO
2355             sums_l(nzb,48,tn) = s1
2356             sums_l(nzb,45,tn) = s2
2357             !$acc end parallel
2358
2359             IF ( cloud_droplets )  THEN
2360
2361                !$OMP DO
2362                !$acc parallel present( pt, q, ql, qsws, rmask, shf, sums_l ) create( s1 )
2363                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2364                DO  i = nxl, nxr
2365                   DO  j =  nys, nyn
2366                      s1 = s1 + ( ( 1.0_wp +                                   &
2367                                    0.61_wp * q(nzb,j,i) - ql(nzb,j,i) ) *     &
2368                                  shf(j,i) + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )
2369                   ENDDO
2370                ENDDO
2371                sums_l(nzb,45,tn) = s1
2372                !$acc end parallel
2373
2374             ENDIF
2375
2376             IF ( cloud_physics )  THEN
2377
2378                !$OMP DO
2379                !$acc parallel present( qsws, rmask, sums_l ) create( s1 )
2380                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2381                DO  i = nxl, nxr
2382                   DO  j =  nys, nyn
2383!
2384!--                   Formula does not work if ql(nzb) /= 0.0
2385                      s1 = s1 + qsws(j,i) * rmask(j,i,sr)   ! w"q" (w"qv")
2386                   ENDDO
2387                ENDDO
2388                sums_l(nzb,51,tn) = s1
2389                !$acc end parallel
2390
2391             ENDIF
2392
2393          ENDIF
2394
2395          IF ( passive_scalar )  THEN
2396
2397             !$OMP DO
2398             !$acc parallel present( qsws, rmask, sums_l ) create( s1 )
2399             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2400             DO  i = nxl, nxr
2401                DO  j =  nys, nyn
2402                   s1 = s1 + qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
2403                ENDDO
2404             ENDDO
2405             sums_l(nzb,48,tn) = s1
2406             !$acc end parallel
2407
2408          ENDIF
2409
2410       ENDIF
2411
2412!
2413!--    Subgridscale fluxes at the top surface
2414       IF ( use_top_fluxes )  THEN
2415
2416          !$OMP DO
2417          !$acc parallel present( rmask, sums_l, tswst, uswst, vswst ) create( s1, s2, s3, s4, s5 )
2418          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5 )
2419          DO  i = nxl, nxr
2420             DO  j =  nys, nyn
2421                s1 = s1 + uswst(j,i) * rmask(j,i,sr)    ! w"u"
2422                s2 = s2 + vswst(j,i) * rmask(j,i,sr)    ! w"v"
2423                s3 = s3 + tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
2424                s4 = s4 + 0.0_wp * rmask(j,i,sr)        ! u"pt"
2425                s5 = s5 + 0.0_wp * rmask(j,i,sr)        ! v"pt"
2426             ENDDO
2427          ENDDO
2428          sums_l(nzt:nzt+1,12,tn) = s1
2429          sums_l(nzt:nzt+1,14,tn) = s2
2430          sums_l(nzt:nzt+1,16,tn) = s3
2431          sums_l(nzt:nzt+1,58,tn) = s4
2432          sums_l(nzt:nzt+1,61,tn) = s5
2433          !$acc end parallel
2434
2435          IF ( ocean )  THEN
2436
2437             !$OMP DO
2438             !$acc parallel present( rmask, saswst, sums_l ) create( s1 )
2439             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2440             DO  i = nxl, nxr
2441                DO  j =  nys, nyn
2442                   s1 = s1 + saswst(j,i) * rmask(j,i,sr)  ! w"sa"
2443                ENDDO
2444             ENDDO
2445             sums_l(nzt,65,tn) = s1
2446             !$acc end parallel
2447
2448          ENDIF
2449
2450          IF ( humidity )  THEN
2451
2452             !$OMP DO
2453             !$acc parallel present( pt, q, qswst, rmask, tswst, sums_l ) create( s1, s2 )
2454             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2455             DO  i = nxl, nxr
2456                DO  j =  nys, nyn
2457                   s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
2458                   s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) * tswst(j,i) +&
2459                                 0.61_wp * pt(nzt,j,i) * qswst(j,i) )
2460                ENDDO
2461             ENDDO
2462             sums_l(nzt,48,tn) = s1
2463             sums_l(nzt,45,tn) = s2
2464             !$acc end parallel
2465
2466             IF ( cloud_droplets )  THEN
2467
2468                !$OMP DO
2469                !$acc parallel present( pt, q, ql, qswst, rmask, tswst, sums_l ) create( s1 )
2470                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2471                DO  i = nxl, nxr
2472                   DO  j =  nys, nyn
2473                      s1 = s1 + ( ( 1.0_wp +                                   &
2474                                    0.61_wp * q(nzt,j,i) - ql(nzt,j,i) ) *     &
2475                                  tswst(j,i) +                                 &
2476                                  0.61_wp * pt(nzt,j,i) * qswst(j,i) )
2477                   ENDDO
2478                ENDDO
2479                sums_l(nzt,45,tn) = s1
2480                !$acc end parallel
2481
2482             ENDIF
2483
2484             IF ( cloud_physics )  THEN
2485
2486                !$OMP DO
2487                !$acc parallel present( qswst, rmask, sums_l ) create( s1 )
2488                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2489                DO  i = nxl, nxr
2490                   DO  j =  nys, nyn
2491!
2492!--                   Formula does not work if ql(nzb) /= 0.0
2493                      s1 = s1 + qswst(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
2494                   ENDDO
2495                ENDDO
2496                sums_l(nzt,51,tn) = s1
2497                !$acc end parallel
2498
2499             ENDIF
2500
2501          ENDIF
2502
2503          IF ( passive_scalar )  THEN
2504
2505             !$OMP DO
2506             !$acc parallel present( qswst, rmask, sums_l ) create( s1 )
2507             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2508             DO  i = nxl, nxr
2509                DO  j =  nys, nyn
2510                   s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
2511                ENDDO
2512             ENDDO
2513             sums_l(nzt,48,tn) = s1
2514             !$acc end parallel
2515
2516          ENDIF
2517
2518       ENDIF
2519
2520!
2521!--    Resolved fluxes (can be computed for all horizontal points)
2522!--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
2523!--    ----  speaking the following k-loop would have to be split up and
2524!--          rearranged according to the staggered grid.
2525       !$acc parallel loop gang present( hom, pt, rflags_invers, rmask, sums_l, u, v, w ) create( s1, s2, s3 )
2526       DO  k = nzb, nzt_diff
2527          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
2528          DO  i = nxl, nxr
2529             DO  j = nys, nyn
2530                ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) + &
2531                                 u(k+1,j,i) - hom(k+1,1,1,sr) )
2532                vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) + &
2533                                 v(k+1,j,i) - hom(k+1,1,2,sr) )
2534                pts = 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) + &
2535                                 pt(k+1,j,i) - hom(k+1,1,4,sr) )
2536!
2537!--             Higher moments
2538                s1 = s1 + pts * w(k,j,i)**2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2539                s2 = s2 + pts**2 * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2540!
2541!--             Energy flux w*e* (has to be adjusted?)
2542                s3 = s3 + w(k,j,i) * 0.5_wp * ( ust**2 + vst**2 + w(k,j,i)**2 )&
2543                                   * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2544             ENDDO
2545          ENDDO
2546          sums_l(k,35,tn) = s1
2547          sums_l(k,36,tn) = s2
2548          sums_l(k,37,tn) = s3
2549       ENDDO
2550       !$acc end parallel loop
2551
2552!
2553!--    Salinity flux and density (density does not belong to here,
2554!--    but so far there is no other suitable place to calculate)
2555       IF ( ocean )  THEN
2556
2557          IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
2558
2559             !$acc parallel loop gang present( hom, rflags_invers, rmask, sa, sums_l, w ) create( s1 )
2560             DO  k = nzb, nzt_diff
2561                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2562                DO  i = nxl, nxr
2563                   DO  j = nys, nyn
2564                      s1 = s1 + 0.5_wp * ( sa(k,j,i)   - hom(k,1,23,sr) +      &
2565                                           sa(k+1,j,i) - hom(k+1,1,23,sr) )    &
2566                                       * w(k,j,i) * rmask(j,i,sr)              & 
2567                                       * rflags_invers(j,i,k+1)
2568                   ENDDO
2569                ENDDO
2570                sums_l(k,66,tn) = s1
2571             ENDDO
2572             !$acc end parallel loop
2573
2574          ENDIF
2575
2576          !$acc parallel loop gang present( rflags_invers, rho, prho, rmask, sums_l ) create( s1, s2 )
2577          DO  k = nzb, nzt_diff
2578             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2579             DO  i = nxl, nxr
2580                DO  j = nys, nyn
2581                   s1 = s1 + rho(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2582                   s2 = s2 + prho(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2583                ENDDO
2584             ENDDO
2585             sums_l(k,64,tn) = s1
2586             sums_l(k,71,tn) = s2
2587          ENDDO
2588          !$acc end parallel loop
2589
2590       ENDIF
2591
2592!
2593!--    Buoyancy flux, water flux, humidity flux, liquid water
2594!--    content, rain drop concentration and rain water content
2595       IF ( humidity )  THEN
2596
2597          IF ( cloud_physics  .OR.  cloud_droplets )  THEN
2598
2599             !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, vpt, w ) create( s1 )
2600             DO  k = nzb, nzt_diff
2601                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2602                DO  i = nxl, nxr
2603                   DO  j = nys, nyn
2604                      s1 = s1 + 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +     &
2605                                           vpt(k+1,j,i) - hom(k+1,1,44,sr) ) * &
2606                                         w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2607                   ENDDO
2608                ENDDO
2609                sums_l(k,46,tn) = s1
2610             ENDDO
2611             !$acc end parallel loop
2612
2613             IF ( .NOT. cloud_droplets )  THEN
2614
2615                !$acc parallel loop gang present( hom, q, ql, rflags_invers, rmask, sums_l, w ) create( s1 )
2616                DO  k = nzb, nzt_diff
2617                   !$acc loop vector collapse( 2 ) reduction( +: s1 )
2618                   DO  i = nxl, nxr
2619                      DO  j = nys, nyn
2620                         s1 = s1 + 0.5_wp * ( ( q(k,j,i)   - ql(k,j,i)   ) - hom(k,1,42,sr) +   &
2621                                              ( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) ) &
2622                                          * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2623                      ENDDO
2624                   ENDDO
2625                   sums_l(k,52,tn) = s1
2626                ENDDO
2627                !$acc end parallel loop
2628
2629                IF ( icloud_scheme == 0  )  THEN
2630
2631                   !$acc parallel loop gang present( qc, ql, rflags_invers, rmask, sums_l ) create( s1, s2 )
2632                   DO  k = nzb, nzt_diff
2633                      !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2634                      DO  i = nxl, nxr
2635                         DO  j = nys, nyn
2636                            s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2637                            s2 = s2 + qc(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2638                         ENDDO
2639                      ENDDO
2640                      sums_l(k,54,tn) = s1
2641                      sums_l(k,75,tn) = s2
2642                   ENDDO
2643                   !$acc end parallel loop
2644
2645                   IF ( precipitation )  THEN
2646
2647                      !$acc parallel loop gang present( nr, qr, prr, rflags_invers, rmask, sums_l ) create( s1, s2, s3 )
2648                      DO  k = nzb, nzt_diff
2649                         !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
2650                         DO  i = nxl, nxr
2651                            DO  j = nys, nyn
2652                               s1 = s1 + nr(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2653                               s2 = s2 + qr(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2654                               s3 = s3 + prr(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2655                            ENDDO
2656                         ENDDO
2657                         sums_l(k,73,tn) = s1
2658                         sums_l(k,74,tn) = s2
2659                         sums_l(k,76,tn) = s3
2660                      ENDDO
2661                      !$acc end parallel loop
2662
2663                   ENDIF
2664
2665                ELSE
2666
2667                   !$acc parallel loop gang present( ql, rflags_invers, rmask, sums_l ) create( s1 )
2668                   DO  k = nzb, nzt_diff
2669                      !$acc loop vector collapse( 2 ) reduction( +: s1 )
2670                      DO  i = nxl, nxr
2671                         DO  j = nys, nyn
2672                            s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2673                         ENDDO
2674                      ENDDO
2675                      sums_l(k,54,tn) = s1
2676                   ENDDO
2677                   !$acc end parallel loop
2678
2679                ENDIF
2680
2681             ELSE
2682
2683                !$acc parallel loop gang present( ql, rflags_invers, rmask, sums_l ) create( s1 )
2684                DO  k = nzb, nzt_diff
2685                   !$acc loop vector collapse( 2 ) reduction( +: s1 )
2686                   DO  i = nxl, nxr
2687                      DO  j = nys, nyn
2688                         s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2689                      ENDDO
2690                   ENDDO
2691                   sums_l(k,54,tn) = s1
2692                ENDDO
2693                !$acc end parallel loop
2694
2695             ENDIF
2696
2697          ELSE
2698
2699             IF( .NOT. ws_scheme_sca  .OR.  sr /= 0 )  THEN
2700
2701                !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, vpt, w ) create( s1 )
2702                DO  k = nzb, nzt_diff
2703                   !$acc loop vector collapse( 2 ) reduction( +: s1 )
2704                   DO  i = nxl, nxr
2705                      DO  j = nys, nyn
2706                         s1 = s1 + 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +   &
2707                                              vpt(k+1,j,i) - hom(k+1,1,44,sr) ) &
2708                                          * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2709                      ENDDO
2710                   ENDDO
2711                   sums_l(k,46,tn) = s1
2712                ENDDO
2713                !$acc end parallel loop
2714
2715             ELSEIF ( ws_scheme_sca  .AND.  sr == 0 )  THEN
2716
2717                !$acc parallel loop present( hom, sums_l )
2718                DO  k = nzb, nzt_diff
2719                   sums_l(k,46,tn) = ( 1.0_wp + 0.61_wp * hom(k,1,41,sr) ) * sums_l(k,17,tn) + &
2720                                                0.61_wp * hom(k,1,4,sr) * sums_l(k,49,tn)
2721                ENDDO
2722                !$acc end parallel loop
2723
2724             ENDIF
2725
2726          ENDIF
2727
2728       ENDIF
2729!
2730!--    Passive scalar flux
2731       IF ( passive_scalar  .AND.  ( .NOT. ws_scheme_sca  .OR.  sr /= 0 ) )  THEN
2732
2733          !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l, w ) create( s1 )
2734          DO  k = nzb, nzt_diff
2735             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2736             DO  i = nxl, nxr
2737                DO  j = nys, nyn
2738                   s1 = s1 + 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +          &
2739                                        q(k+1,j,i) - hom(k+1,1,41,sr) )        &
2740                                    * w(k,j,i) * rmask(j,i,sr)                 &
2741                                    * rflags_invers(j,i,k+1)
2742                ENDDO
2743             ENDDO
2744             sums_l(k,49,tn) = s1
2745          ENDDO
2746          !$acc end parallel loop
2747
2748       ENDIF
2749
2750!
2751!--    For speed optimization fluxes which have been computed in part directly
2752!--    inside the WS advection routines are treated seperatly
2753!--    Momentum fluxes first:
2754       IF ( .NOT. ws_scheme_mom  .OR.  sr /= 0  )  THEN
2755
2756          !$OMP DO
2757          !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, u, v, w ) create( s1, s2 )
2758          DO  k = nzb, nzt_diff
2759             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2760             DO  i = nxl, nxr
2761                DO  j = nys, nyn
2762                   ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +               &
2763                                    u(k+1,j,i) - hom(k+1,1,1,sr) )
2764                   vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +               &
2765                                    v(k+1,j,i) - hom(k+1,1,2,sr) )
2766!
2767!--                Momentum flux w*u*
2768                   s1 = s1 + 0.5_wp * ( w(k,j,i-1) + w(k,j,i) )                &
2769                                    * ust * rmask(j,i,sr)                      &
2770                                    * rflags_invers(j,i,k+1)
2771!
2772!--                Momentum flux w*v*
2773                   s2 = s2 + 0.5_wp * ( w(k,j-1,i) + w(k,j,i) )                &
2774                                    * vst * rmask(j,i,sr)                      & 
2775                                    * rflags_invers(j,i,k+1)
2776                ENDDO
2777             ENDDO
2778             sums_l(k,13,tn) = s1
2779             sums_l(k,15,tn) = s1
2780          ENDDO
2781          !$acc end parallel loop
2782
2783       ENDIF
2784
2785       IF ( .NOT. ws_scheme_sca  .OR.  sr /= 0 )  THEN
2786
2787          !$OMP DO
2788          !$acc parallel loop gang present( hom, pt, rflags_invers, rmask, sums_l, w ) create( s1 )
2789          DO  k = nzb, nzt_diff
2790             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2791             DO  i = nxl, nxr
2792                DO  j = nys, nyn
2793!
2794!--                Vertical heat flux
2795                   s1 = s1 + 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) +          &
2796                                        pt(k+1,j,i) - hom(k+1,1,4,sr) )        &
2797                                    * w(k,j,i) * rmask(j,i,sr)                 & 
2798                                    * rflags_invers(j,i,k+1)
2799                ENDDO
2800             ENDDO
2801             sums_l(k,17,tn) = s1
2802          ENDDO
2803          !$acc end parallel loop
2804
2805          IF ( humidity )  THEN
2806
2807             !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l, w ) create( s1 )
2808             DO  k = nzb, nzt_diff
2809                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2810                DO  i = nxl, nxr
2811                   DO  j = nys, nyn
2812                      s1 = s1 + 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +       &
2813                                           q(k+1,j,i) - hom(k+1,1,41,sr) )     &
2814                                       * w(k,j,i) * rmask(j,i,sr)              &
2815                                       * rflags_invers(j,i,k+1)
2816                   ENDDO
2817                ENDDO
2818                sums_l(k,49,tn) = s1
2819             ENDDO
2820             !$acc end parallel loop
2821
2822          ENDIF
2823
2824       ENDIF
2825
2826
2827!
2828!--    Density at top follows Neumann condition
2829       IF ( ocean )  THEN
2830          !$acc parallel present( sums_l )
2831          sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
2832          sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
2833          !$acc end parallel
2834       ENDIF
2835
2836!
2837!--    Divergence of vertical flux of resolved scale energy and pressure
2838!--    fluctuations as well as flux of pressure fluctuation itself (68).
2839!--    First calculate the products, then the divergence.
2840!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
2841       IF ( hom(nzb+1,2,55,0) /= 0.0_wp  .OR.  hom(nzb+1,2,68,0) /= 0.0_wp )  THEN
2842
2843          STOP '+++ openACC porting for vertical flux div of resolved scale TKE in flow_statistics is still missing'
2844          sums_ll = 0.0_wp  ! local array
2845
2846          !$OMP DO
2847          DO  i = nxl, nxr
2848             DO  j = nys, nyn
2849                DO  k = nzb_s_inner(j,i)+1, nzt
2850
2851                   sums_ll(k,1) = sums_ll(k,1) + 0.5_wp * w(k,j,i) * (         &
2852                  ( 0.25_wp * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1)    &
2853                              - 0.5_wp * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) )   &
2854                              ) )**2                                           &
2855                + ( 0.25_wp * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i)    &
2856                              - 0.5_wp * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) )   &
2857                              ) )**2                                           &
2858                + w(k,j,i)**2                                        )
2859
2860                   sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i)             &
2861                                               * ( p(k,j,i) + p(k+1,j,i) )
2862
2863                ENDDO
2864             ENDDO
2865          ENDDO
2866          sums_ll(0,1)     = 0.0_wp    ! because w is zero at the bottom
2867          sums_ll(nzt+1,1) = 0.0_wp
2868          sums_ll(0,2)     = 0.0_wp
2869          sums_ll(nzt+1,2) = 0.0_wp
2870
2871          DO  k = nzb+1, nzt
2872             sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
2873             sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
2874             sums_l(k,68,tn) = sums_ll(k,2)
2875          ENDDO
2876          sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
2877          sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
2878          sums_l(nzb,68,tn) = 0.0_wp    ! because w* = 0 at nzb
2879
2880       ENDIF
2881
2882!
2883!--    Divergence of vertical flux of SGS TKE and the flux itself (69)
2884       IF ( hom(nzb+1,2,57,0) /= 0.0_wp  .OR.  hom(nzb+1,2,69,0) /= 0.0_wp )  THEN
2885
2886          STOP '+++ openACC porting for vertical flux div of SGS TKE in flow_statistics is still missing'
2887          !$OMP DO
2888          DO  i = nxl, nxr
2889             DO  j = nys, nyn
2890                DO  k = nzb_s_inner(j,i)+1, nzt
2891
2892                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5_wp * (              &
2893                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
2894                 - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
2895                                                                ) * ddzw(k)
2896
2897                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5_wp * (              &
2898                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
2899                                                                )
2900
2901                ENDDO
2902             ENDDO
2903          ENDDO
2904          sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
2905          sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
2906
2907       ENDIF
2908
2909!
2910!--    Horizontal heat fluxes (subgrid, resolved, total).
2911!--    Do it only, if profiles shall be plotted.
2912       IF ( hom(nzb+1,2,58,0) /= 0.0_wp ) THEN
2913
2914          STOP '+++ openACC porting for horizontal flux calculation in flow_statistics is still missing'
2915          !$OMP DO
2916          DO  i = nxl, nxr
2917             DO  j = nys, nyn
2918                DO  k = nzb_s_inner(j,i)+1, nzt
2919!
2920!--                Subgrid horizontal heat fluxes u"pt", v"pt"
2921                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5_wp *                &
2922                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
2923                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
2924                                                 * ddx * rmask(j,i,sr)
2925                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
2926                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
2927                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
2928                                                 * ddy * rmask(j,i,sr)
2929!
2930!--                Resolved horizontal heat fluxes u*pt*, v*pt*
2931                   sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
2932                                     ( u(k,j,i) - hom(k,1,1,sr) ) * 0.5_wp *   &
2933                                     ( pt(k,j,i-1) - hom(k,1,4,sr) +           &
2934                                       pt(k,j,i)   - hom(k,1,4,sr) )
2935                   pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
2936                                    pt(k,j,i)   - hom(k,1,4,sr) )
2937                   sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
2938                                     ( v(k,j,i) - hom(k,1,2,sr) ) * 0.5_wp *   &
2939                                     ( pt(k,j-1,i) - hom(k,1,4,sr) +           &
2940                                       pt(k,j,i)   - hom(k,1,4,sr) )
2941                ENDDO
2942             ENDDO
2943          ENDDO
2944!
2945!--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
2946          sums_l(nzb,58,tn) = 0.0_wp
2947          sums_l(nzb,59,tn) = 0.0_wp
2948          sums_l(nzb,60,tn) = 0.0_wp
2949          sums_l(nzb,61,tn) = 0.0_wp
2950          sums_l(nzb,62,tn) = 0.0_wp
2951          sums_l(nzb,63,tn) = 0.0_wp
2952
2953       ENDIF
2954
2955!
2956!--    Collect current large scale advection and subsidence tendencies for
2957!--    data output
2958       IF ( large_scale_forcing  .AND.  ( simulated_time .GT. 0.0_wp ) )  THEN
2959!
2960!--       Interpolation in time of LSF_DATA
2961          nt = 1
2962          DO WHILE ( simulated_time - dt_3d > time_vert(nt) )
2963             nt = nt + 1
2964          ENDDO
2965          IF ( simulated_time - dt_3d /= time_vert(nt) )  THEN
2966            nt = nt - 1
2967          ENDIF
2968
2969          fac = ( simulated_time - dt_3d - time_vert(nt) )                     &
2970                / ( time_vert(nt+1)-time_vert(nt) )
2971
2972
2973          DO  k = nzb, nzt
2974             sums_ls_l(k,0) = td_lsa_lpt(k,nt)                                 &
2975                              + fac * ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )
2976             sums_ls_l(k,1) = td_lsa_q(k,nt)                                   &
2977                              + fac * ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) )
2978          ENDDO
2979
2980          sums_ls_l(nzt+1,0) = sums_ls_l(nzt,0)
2981          sums_ls_l(nzt+1,1) = sums_ls_l(nzt,1)
2982
2983          IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
2984
2985             DO  k = nzb, nzt
2986                sums_ls_l(k,2) = td_sub_lpt(k,nt) + fac *                      &
2987                                 ( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )
2988                sums_ls_l(k,3) = td_sub_q(k,nt) + fac *                        &
2989                                 ( td_sub_q(k,nt+1) - td_sub_q(k,nt) )
2990             ENDDO
2991
2992             sums_ls_l(nzt+1,2) = sums_ls_l(nzt,2)
2993             sums_ls_l(nzt+1,3) = sums_ls_l(nzt,3)
2994
2995          ENDIF
2996
2997       ENDIF
2998
2999
3000       IF ( land_surface )  THEN
3001          !$OMP DO
3002          DO  i = nxl, nxr
3003             DO  j =  nys, nyn
3004                DO  k = nzb_soil, nzt_soil
3005                   sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil(k,j,i) * rmask(j,i,sr)
3006                   sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil(k,j,i) * rmask(j,i,sr)
3007                ENDDO
3008             ENDDO
3009          ENDDO
3010       ENDIF
3011
3012
3013       IF ( radiation .AND. radiation_scheme == 'rrtmg' )  THEN
3014          !$OMP DO
3015          DO  i = nxl, nxr
3016             DO  j =  nys, nyn
3017                DO  k = nzb_s_inner(j,i)+1, nzt+1
3018                   sums_l(k,102,tn)  = sums_l(k,102,tn)  + rad_lw_in(k,j,i) * rmask(j,i,sr)
3019                   sums_l(k,103,tn)  = sums_l(k,103,tn)  + rad_lw_out(k,j,i) * rmask(j,i,sr)
3020                   sums_l(k,104,tn)  = sums_l(k,104,tn)  + rad_sw_in(k,j,i) * rmask(j,i,sr)
3021                   sums_l(k,105,tn)  = sums_l(k,105,tn)  + rad_sw_out(k,j,i) * rmask(j,i,sr)
3022                ENDDO
3023             ENDDO
3024          ENDDO
3025       ENDIF
3026
3027!
3028!--    Calculate the user-defined profiles
3029       CALL user_statistics( 'profiles', sr, tn )
3030       !$OMP END PARALLEL
3031
3032!
3033!--    Summation of thread sums
3034       IF ( threads_per_task > 1 )  THEN
3035          STOP '+++ openACC porting for threads_per_task > 1 in flow_statistics is still missing'
3036          DO  i = 1, threads_per_task-1
3037             sums_l(:,3,0)          = sums_l(:,3,0) + sums_l(:,3,i)
3038             sums_l(:,4:40,0)       = sums_l(:,4:40,0) + sums_l(:,4:40,i)
3039             sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
3040                                      sums_l(:,45:pr_palm,i)
3041             IF ( max_pr_user > 0 )  THEN
3042                sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
3043                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
3044                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
3045             ENDIF
3046          ENDDO
3047       ENDIF
3048
3049       !$acc update host( hom, sums, sums_l )
3050
3051#if defined( __parallel )
3052
3053!
3054!--    Compute total sum from local sums
3055       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
3056       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, &
3057                           MPI_SUM, comm2d, ierr )
3058       IF ( large_scale_forcing )  THEN
3059          CALL MPI_ALLREDUCE( sums_ls_l(nzb,2), sums(nzb,83), ngp_sums_ls,     &
3060                              MPI_REAL, MPI_SUM, comm2d, ierr )
3061       ENDIF
3062#else
3063       sums = sums_l(:,:,0)
3064       IF ( large_scale_forcing )  THEN
3065          sums(:,81:88) = sums_ls_l
3066       ENDIF
3067#endif
3068
3069!
3070!--    Final values are obtained by division by the total number of grid points
3071!--    used for summation. After that store profiles.
3072!--    Profiles:
3073       DO  k = nzb, nzt+1
3074          sums(k,3)               = sums(k,3)           / ngp_2dh(sr)
3075          sums(k,8:11)            = sums(k,8:11)        / ngp_2dh_s_inner(k,sr)
3076          sums(k,12:22)           = sums(k,12:22)       / ngp_2dh(sr)
3077          sums(k,23:29)           = sums(k,23:29)       / ngp_2dh_s_inner(k,sr)
3078          sums(k,30:32)           = sums(k,30:32)       / ngp_2dh(sr)
3079          sums(k,33:34)           = sums(k,33:34)       / ngp_2dh_s_inner(k,sr)
3080          sums(k,35:39)           = sums(k,35:39)       / ngp_2dh(sr)
3081          sums(k,40)              = sums(k,40)          / ngp_2dh_s_inner(k,sr)
3082          sums(k,45:53)           = sums(k,45:53)       / ngp_2dh(sr)
3083          sums(k,54)              = sums(k,54)          / ngp_2dh_s_inner(k,sr)
3084          sums(k,55:63)           = sums(k,55:63)       / ngp_2dh(sr)
3085          sums(k,64)              = sums(k,64)          / ngp_2dh_s_inner(k,sr)
3086          sums(k,70:80)           = sums(k,70:80)       / ngp_2dh_s_inner(k,sr)
3087          sums(k,81:88)           = sums(k,81:88)       / ngp_2dh(sr)
3088          sums(k,89:105)           = sums(k,89:105)     / ngp_2dh(sr)
3089          sums(k,106:pr_palm-2)    = sums(k,106:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
3090       ENDDO
3091
3092!--    Upstream-parts
3093       sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
3094!--    u* and so on
3095!--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
3096!--    size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
3097!--    above the topography, they are being divided by ngp_2dh(sr)
3098       sums(nzb:nzb+3,pr_palm)    = sums(nzb:nzb+3,pr_palm)    / &
3099                                    ngp_2dh(sr)
3100       sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)       / &    ! qs
3101                                    ngp_2dh(sr)
3102!--    eges, e*
3103       sums(nzb+4:nzb+5,pr_palm)  = sums(nzb+4:nzb+5,pr_palm)  / &
3104                                    ngp_3d(sr)
3105!--    Old and new divergence
3106       sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
3107                                    ngp_3d_inner(sr)
3108
3109!--    User-defined profiles
3110       IF ( max_pr_user > 0 )  THEN
3111          DO  k = nzb, nzt+1
3112             sums(k,pr_palm+1:pr_palm+max_pr_user) = &
3113                                    sums(k,pr_palm+1:pr_palm+max_pr_user) / &
3114                                    ngp_2dh_s_inner(k,sr)
3115          ENDDO
3116       ENDIF
3117
3118!
3119!--    Collect horizontal average in hom.
3120!--    Compute deduced averages (e.g. total heat flux)
3121       hom(:,1,3,sr)  = sums(:,3)      ! w
3122       hom(:,1,8,sr)  = sums(:,8)      ! e     profiles 5-7 are initial profiles
3123       hom(:,1,9,sr)  = sums(:,9)      ! km
3124       hom(:,1,10,sr) = sums(:,10)     ! kh
3125       hom(:,1,11,sr) = sums(:,11)     ! l
3126       hom(:,1,12,sr) = sums(:,12)     ! w"u"
3127       hom(:,1,13,sr) = sums(:,13)     ! w*u*
3128       hom(:,1,14,sr) = sums(:,14)     ! w"v"
3129       hom(:,1,15,sr) = sums(:,15)     ! w*v*
3130       hom(:,1,16,sr) = sums(:,16)     ! w"pt"
3131       hom(:,1,17,sr) = sums(:,17)     ! w*pt*
3132       hom(:,1,18,sr) = sums(:,16) + sums(:,17)    ! wpt
3133       hom(:,1,19,sr) = sums(:,12) + sums(:,13)    ! wu
3134       hom(:,1,20,sr) = sums(:,14) + sums(:,15)    ! wv
3135       hom(:,1,21,sr) = sums(:,21)     ! w*pt*BC
3136       hom(:,1,22,sr) = sums(:,16) + sums(:,21)    ! wptBC
3137                                       ! profile 24 is initial profile (sa)
3138                                       ! profiles 25-29 left empty for initial
3139                                       ! profiles
3140       hom(:,1,30,sr) = sums(:,30)     ! u*2
3141       hom(:,1,31,sr) = sums(:,31)     ! v*2
3142       hom(:,1,32,sr) = sums(:,32)     ! w*2
3143       hom(:,1,33,sr) = sums(:,33)     ! pt*2
3144       hom(:,1,34,sr) = sums(:,34)     ! e*
3145       hom(:,1,35,sr) = sums(:,35)     ! w*2pt*
3146       hom(:,1,36,sr) = sums(:,36)     ! w*pt*2
3147       hom(:,1,37,sr) = sums(:,37)     ! w*e*
3148       hom(:,1,38,sr) = sums(:,38)     ! w*3
3149       hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20_wp )**1.5_wp   ! Sw
3150       hom(:,1,40,sr) = sums(:,40)     ! p
3151       hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
3152       hom(:,1,46,sr) = sums(:,46)     ! w*vpt*
3153       hom(:,1,47,sr) = sums(:,45) + sums(:,46)    ! wvpt
3154       hom(:,1,48,sr) = sums(:,48)     ! w"q" (w"qv")
3155       hom(:,1,49,sr) = sums(:,49)     ! w*q* (w*qv*)
3156       hom(:,1,50,sr) = sums(:,48) + sums(:,49)    ! wq (wqv)
3157       hom(:,1,51,sr) = sums(:,51)     ! w"qv"
3158       hom(:,1,52,sr) = sums(:,52)     ! w*qv*
3159       hom(:,1,53,sr) = sums(:,52) + sums(:,51)    ! wq (wqv)
3160       hom(:,1,54,sr) = sums(:,54)     ! ql
3161       hom(:,1,55,sr) = sums(:,55)     ! w*u*u*/dz
3162       hom(:,1,56,sr) = sums(:,56)     ! w*p*/dz
3163       hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho )/dz
3164       hom(:,1,58,sr) = sums(:,58)     ! u"pt"
3165       hom(:,1,59,sr) = sums(:,59)     ! u*pt*
3166       hom(:,1,60,sr) = sums(:,58) + sums(:,59)    ! upt_t
3167       hom(:,1,61,sr) = sums(:,61)     ! v"pt"
3168       hom(:,1,62,sr) = sums(:,62)     ! v*pt*
3169       hom(:,1,63,sr) = sums(:,61) + sums(:,62)    ! vpt_t
3170       hom(:,1,64,sr) = sums(:,64)     ! rho
3171       hom(:,1,65,sr) = sums(:,65)     ! w"sa"
3172       hom(:,1,66,sr) = sums(:,66)     ! w*sa*
3173       hom(:,1,67,sr) = sums(:,65) + sums(:,66)    ! wsa
3174       hom(:,1,68,sr) = sums(:,68)     ! w*p*
3175       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho
3176       hom(:,1,70,sr) = sums(:,70)     ! q*2
3177       hom(:,1,71,sr) = sums(:,71)     ! prho
3178       hom(:,1,72,sr) = hyp * 1E-4_wp     ! hyp in dbar
3179       hom(:,1,73,sr) = sums(:,73)     ! nr
3180       hom(:,1,74,sr) = sums(:,74)     ! qr
3181       hom(:,1,75,sr) = sums(:,75)     ! qc
3182       hom(:,1,76,sr) = sums(:,76)     ! prr (precipitation rate)
3183                                       ! 77 is initial density profile
3184       hom(:,1,78,sr) = ug             ! ug
3185       hom(:,1,79,sr) = vg             ! vg
3186       hom(:,1,80,sr) = w_subs         ! w_subs
3187
3188       IF ( large_scale_forcing )  THEN
3189          hom(:,1,81,sr) = sums_ls_l(:,0)          ! td_lsa_lpt
3190          hom(:,1,82,sr) = sums_ls_l(:,1)          ! td_lsa_q
3191          IF ( use_subsidence_tendencies )  THEN
3192             hom(:,1,83,sr) = sums_ls_l(:,2)       ! td_sub_lpt
3193             hom(:,1,84,sr) = sums_ls_l(:,3)       ! td_sub_q
3194          ELSE
3195             hom(:,1,83,sr) = sums(:,83)           ! td_sub_lpt
3196             hom(:,1,84,sr) = sums(:,84)           ! td_sub_q
3197          ENDIF
3198          hom(:,1,85,sr) = sums(:,85)              ! td_nud_lpt
3199          hom(:,1,86,sr) = sums(:,86)              ! td_nud_q
3200          hom(:,1,87,sr) = sums(:,87)              ! td_nud_u
3201          hom(:,1,88,sr) = sums(:,88)              ! td_nud_v
3202       END IF
3203
3204       hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
3205                                       ! upstream-parts u_x, u_y, u_z, v_x,
3206                                       ! v_y, usw. (in last but one profile)
3207       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
3208                                       ! u*, w'u', w'v', t* (in last profile)
3209
3210       IF ( max_pr_user > 0 )  THEN    ! user-defined profiles
3211          hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
3212                               sums(:,pr_palm+1:pr_palm+max_pr_user)
3213       ENDIF
3214
3215!
3216!--    Determine the boundary layer height using two different schemes.
3217!--    First scheme: Starting from the Earth's (Ocean's) surface, look for the
3218!--    first relative minimum (maximum) of the total heat flux.
3219!--    The corresponding height is assumed as the boundary layer height, if it
3220!--    is less than 1.5 times the height where the heat flux becomes negative
3221!--    (positive) for the first time.
3222       z_i(1) = 0.0_wp
3223       first = .TRUE.
3224
3225       IF ( ocean )  THEN
3226          DO  k = nzt, nzb+1, -1
3227             IF ( first .AND. hom(k,1,18,sr) < 0.0_wp                          &
3228                .AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp )  THEN
3229                first = .FALSE.
3230                height = zw(k)
3231             ENDIF
3232             IF ( hom(k,1,18,sr) < 0.0_wp  .AND.                               &
3233                  abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND.                        &
3234                  hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
3235                IF ( zw(k) < 1.5_wp * height )  THEN
3236                   z_i(1) = zw(k)
3237                ELSE
3238                   z_i(1) = height
3239                ENDIF
3240                EXIT
3241             ENDIF
3242          ENDDO
3243       ELSE
3244          DO  k = nzb, nzt-1
3245             IF ( first .AND. hom(k,1,18,sr) < 0.0_wp                          &
3246                .AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp )  THEN
3247                first = .FALSE.
3248                height = zw(k)
3249             ENDIF
3250             IF ( hom(k,1,18,sr) < 0.0  .AND. &
3251                  abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND. &
3252                  hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
3253                IF ( zw(k) < 1.5_wp * height )  THEN
3254                   z_i(1) = zw(k)
3255                ELSE
3256                   z_i(1) = height
3257                ENDIF
3258                EXIT
3259             ENDIF
3260          ENDDO
3261       ENDIF
3262
3263!
3264!--    Second scheme: Gradient scheme from Sullivan et al. (1998), modified
3265!--    by Uhlenbrock(2006). The boundary layer height is the height with the
3266!--    maximal local temperature gradient: starting from the second (the last
3267!--    but one) vertical gridpoint, the local gradient must be at least
3268!--    0.2K/100m and greater than the next four gradients.
3269!--    WARNING: The threshold value of 0.2K/100m must be adjusted for the
3270!--             ocean case!
3271       z_i(2) = 0.0_wp
3272       DO  k = nzb+1, nzt+1
3273          dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
3274       ENDDO
3275       dptdz_threshold = 0.2_wp / 100.0_wp
3276
3277       IF ( ocean )  THEN
3278          DO  k = nzt+1, nzb+5, -1
3279             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
3280                  dptdz(k) > dptdz(k-1)  .AND.  dptdz(k) > dptdz(k-2)  .AND.  &
3281                  dptdz(k) > dptdz(k-3)  .AND.  dptdz(k) > dptdz(k-4) )  THEN
3282                z_i(2) = zw(k-1)
3283                EXIT
3284             ENDIF
3285          ENDDO
3286       ELSE
3287          DO  k = nzb+1, nzt-3
3288             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
3289                  dptdz(k) > dptdz(k+1)  .AND.  dptdz(k) > dptdz(k+2)  .AND.  &
3290                  dptdz(k) > dptdz(k+3)  .AND.  dptdz(k) > dptdz(k+4) )  THEN
3291                z_i(2) = zw(k-1)
3292                EXIT
3293             ENDIF
3294          ENDDO
3295       ENDIF
3296
3297       hom(nzb+6,1,pr_palm,sr) = z_i(1)
3298       hom(nzb+7,1,pr_palm,sr) = z_i(2)
3299
3300!
3301!--    Computation of both the characteristic vertical velocity and
3302!--    the characteristic convective boundary layer temperature.
3303!--    The horizontal average at nzb+1 is input for the average temperature.
3304       IF ( hom(nzb,1,18,sr) > 0.0_wp .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8_wp  &
3305           .AND.  z_i(1) /= 0.0_wp )  THEN
3306          hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) *                 &
3307                                       hom(nzb,1,18,sr) *                      &
3308                                       ABS( z_i(1) ) )**0.333333333_wp
3309!--       so far this only works if Prandtl layer is used
3310          hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
3311       ELSE
3312          hom(nzb+8,1,pr_palm,sr)  = 0.0_wp
3313          hom(nzb+11,1,pr_palm,sr) = 0.0_wp
3314       ENDIF
3315
3316!
3317!--    Collect the time series quantities
3318       ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)     ! E
3319       ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)     ! E*
3320       ts_value(3,sr) = dt_3d
3321       ts_value(4,sr) = hom(nzb,1,pr_palm,sr)       ! u*
3322       ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr)     ! th*
3323       ts_value(6,sr) = u_max
3324       ts_value(7,sr) = v_max
3325       ts_value(8,sr) = w_max
3326       ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr)    ! new divergence
3327       ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr)    ! old Divergence
3328       ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr)    ! z_i(1)
3329       ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr)    ! z_i(2)
3330       ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr)    ! w*
3331       ts_value(14,sr) = hom(nzb,1,16,sr)           ! w'pt'   at k=0
3332       ts_value(15,sr) = hom(nzb+1,1,16,sr)         ! w'pt'   at k=1
3333       ts_value(16,sr) = hom(nzb+1,1,18,sr)         ! wpt     at k=1
3334       ts_value(17,sr) = hom(nzb,1,4,sr)            ! pt(0)
3335       ts_value(18,sr) = hom(nzb+1,1,4,sr)          ! pt(zp)
3336       ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr)    ! u'w'    at k=0
3337       ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)    ! v'w'    at k=0
3338       ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
3339
3340       IF ( ts_value(5,sr) /= 0.0_wp )  THEN
3341          ts_value(22,sr) = ts_value(4,sr)**2_wp / &
3342                            ( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
3343       ELSE
3344          ts_value(22,sr) = 10000.0_wp
3345       ENDIF
3346
3347       ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
3348
3349!
3350!--    Collect land surface model timeseries
3351       IF ( land_surface )  THEN
3352          ts_value(dots_soil  ,sr) = hom(nzb,1,93,sr)           ! ghf_eb
3353          ts_value(dots_soil+1,sr) = hom(nzb,1,94,sr)           ! shf_eb
3354          ts_value(dots_soil+2,sr) = hom(nzb,1,95,sr)           ! qsws_eb
3355          ts_value(dots_soil+3,sr) = hom(nzb,1,96,sr)           ! qsws_liq_eb
3356          ts_value(dots_soil+4,sr) = hom(nzb,1,97,sr)           ! qsws_soil_eb
3357          ts_value(dots_soil+5,sr) = hom(nzb,1,98,sr)           ! qsws_veg_eb
3358          ts_value(dots_soil+6,sr) = hom(nzb,1,99,sr)           ! r_a
3359          ts_value(dots_soil+7,sr) = hom(nzb,1,100,sr)          ! r_s
3360       ENDIF
3361!
3362!--    Collect radiation model timeseries
3363       IF ( radiation )  THEN
3364          ts_value(dots_rad,sr)   = hom(nzb,1,101,sr)          ! rad_net
3365          ts_value(dots_rad+1,sr) = hom(nzb,1,102,sr)          ! rad_lw_in
3366          ts_value(dots_rad+2,sr) = hom(nzb,1,103,sr)          ! rad_lw_out
3367          ts_value(dots_rad+3,sr) = hom(nzb,1,104,sr)          ! rad_lw_in
3368          ts_value(dots_rad+4,sr) = hom(nzb,1,105,sr)          ! rad_lw_out
3369
3370#if defined ( __rrtmg )
3371          IF ( radiation_scheme == 'rrtmg' )  THEN
3372             ts_value(dots_rad+5,sr) = hom(nzb,1,106,sr)          ! rrtm_aldif
3373             ts_value(dots_rad+6,sr) = hom(nzb,1,107,sr)          ! rrtm_aldir
3374             ts_value(dots_rad+7,sr) = hom(nzb,1,108,sr)          ! rrtm_asdif
3375             ts_value(dots_rad+8,sr) = hom(nzb,1,109,sr)          ! rrtm_asdir
3376          ENDIF
3377#endif
3378
3379       ENDIF
3380
3381!
3382!--    Calculate additional statistics provided by the user interface
3383       CALL user_statistics( 'time_series', sr, 0 )
3384
3385    ENDDO    ! loop of the subregions
3386
3387    !$acc end data
3388
3389!
3390!-- If required, sum up horizontal averages for subsequent time averaging
3391    IF ( do_sum )  THEN
3392       IF ( average_count_pr == 0 )  hom_sum = 0.0_wp
3393       hom_sum = hom_sum + hom(:,1,:,:)
3394       average_count_pr = average_count_pr + 1
3395       do_sum = .FALSE.
3396    ENDIF
3397
3398!
3399!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
3400!-- This flag is reset after each time step in time_integration.
3401    flow_statistics_called = .TRUE.
3402
3403    CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
3404
3405
3406 END SUBROUTINE flow_statistics
3407#endif
Note: See TracBrowser for help on using the repository browser.