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

Last change on this file since 1557 was 1557, checked in by suehring, 6 years ago

Enable monotone advection for scalars in combination with fifth-order scheme using monotonic limiter.

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