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

Last change on this file since 1985 was 1977, checked in by maronga, 8 years ago

last commit documented

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