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

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

minor bugfixes for radiation model. bugfix in subjob

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