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

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

cpp-switches removed + cpp-bugfixes + zero-settings for velocities inside topography re-activated

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