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

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

last commit documented

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