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

Last change on this file since 1682 was 1682, checked in by knoop, 6 years ago

Code annotations made doxygen readable

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