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

Last change on this file since 1784 was 1784, checked in by raasch, 5 years ago

last commit documented

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