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

Last change on this file since 1572 was 1572, checked in by maronga, 9 years ago

last commit documented

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