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

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

last commit documented

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