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

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

last commit documented

  • Property svn:keywords set to Id
File size: 138.8 KB
Line 
1#if ! defined( __openacc )
2 SUBROUTINE flow_statistics
3
4!--------------------------------------------------------------------------------!
5! This file is part of PALM.
6!
7! PALM is free software: you can redistribute it and/or modify it under the terms
8! of the GNU General Public License as published by the Free Software Foundation,
9! either version 3 of the License, or (at your option) any later version.
10!
11! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
12! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
13! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
14!
15! You should have received a copy of the GNU General Public License along with
16! PALM. If not, see <http://www.gnu.org/licenses/>.
17!
18! Copyright 1997-2014 Leibniz Universitaet Hannover
19!--------------------------------------------------------------------------------!
20!
21! Current revisions:
22! -----------------
23!
24!
25! Former revisions:
26! -----------------
27! $Id: flow_statistics.f90 1558 2015-03-05 16:51:41Z suehring $
28!
29! 1557 2015-03-05 16:43:04Z suehring
30! Adjustments for monotonic limiter
31!
32! 1555 2015-03-04 17:44:27Z maronga
33! Added output of r_a and r_s.
34!
35! 1551 2015-03-03 14:18:16Z maronga
36! Added suppport for land surface model and radiation model output.
37!
38! 1498 2014-12-03 14:09:51Z suehring
39! Comments added
40!
41! 1482 2014-10-18 12:34:45Z raasch
42! missing ngp_sums_ls added in accelerator version
43!
44! 1450 2014-08-21 07:31:51Z heinze
45! bugfix: calculate fac only for simulated_time >= 0.0
46!
47! 1396 2014-05-06 13:37:41Z raasch
48! bugfix: "copyin" replaced by "update device" in openacc-branch
49!
50! 1386 2014-05-05 13:55:30Z boeske
51! bugfix: simulated time before the last timestep is needed for the correct
52! calculation of the profiles of large scale forcing tendencies
53!
54! 1382 2014-04-30 12:15:41Z boeske
55! Renamed variables which store large scale forcing tendencies
56! pt_lsa -> td_lsa_lpt, pt_subs -> td_sub_lpt,
57! q_lsa  -> td_lsa_q,   q_subs  -> td_sub_q,
58! added Neumann boundary conditions for profile data output of large scale
59! advection and subsidence terms at nzt+1
60!
61! 1374 2014-04-25 12:55:07Z raasch
62! bugfix: syntax errors removed from openacc-branch
63! missing variables added to ONLY-lists
64!
65! 1365 2014-04-22 15:03:56Z boeske
66! Output of large scale advection, large scale subsidence and nudging tendencies
67! +sums_ls_l, ngp_sums_ls, use_subsidence_tendencies
68!
69! 1353 2014-04-08 15:21:23Z heinze
70! REAL constants provided with KIND-attribute
71!
72! 1322 2014-03-20 16:38:49Z raasch
73! REAL constants defined as wp-kind
74!
75! 1320 2014-03-20 08:40:49Z raasch
76! ONLY-attribute added to USE-statements,
77! kind-parameters added to all INTEGER and REAL declaration statements,
78! kinds are defined in new module kinds,
79! revision history before 2012 removed,
80! comment fields (!:) to be used for variable explanations added to
81! all variable declaration statements
82!
83! 1299 2014-03-06 13:15:21Z heinze
84! Output of large scale vertical velocity w_subs
85!
86! 1257 2013-11-08 15:18:40Z raasch
87! openacc "end parallel" replaced by "end parallel loop"
88!
89! 1241 2013-10-30 11:36:58Z heinze
90! Output of ug and vg
91!
92! 1221 2013-09-10 08:59:13Z raasch
93! ported for openACC in separate #else branch
94!
95! 1179 2013-06-14 05:57:58Z raasch
96! comment for profile 77 added
97!
98! 1115 2013-03-26 18:16:16Z hoffmann
99! ql is calculated by calc_liquid_water_content
100!
101! 1111 2013-03-08 23:54:10Z raasch
102! openACC directive added
103!
104! 1053 2012-11-13 17:11:03Z hoffmann
105! additions for two-moment cloud physics scheme:
106! +nr, qr, qc, prr
107!
108! 1036 2012-10-22 13:43:42Z raasch
109! code put under GPL (PALM 3.9)
110!
111! 1007 2012-09-19 14:30:36Z franke
112! Calculation of buoyancy flux for humidity in case of WS-scheme is now using
113! turbulent fluxes of WS-scheme
114! Bugfix: Calculation of subgridscale buoyancy flux for humidity and cloud
115! droplets at nzb and nzt added
116!
117! 801 2012-01-10 17:30:36Z suehring
118! Calculation of turbulent fluxes in advec_ws is now thread-safe.
119!
120! Revision 1.1  1997/08/11 06:15:17  raasch
121! Initial revision
122!
123!
124! Description:
125! ------------
126! Compute average profiles and further average flow quantities for the different
127! user-defined (sub-)regions. The region indexed 0 is the total model domain.
128!
129! NOTE: For simplicity, nzb_s_inner and nzb_diff_s_inner are being used as a
130! ----  lower vertical index for k-loops for all variables, although strictly
131! speaking the k-loops would have to be split up according to the staggered
132! grid. However, this implies no error since staggered velocity components are
133! zero at the walls and inside buildings.
134!------------------------------------------------------------------------------!
135
136    USE arrays_3d,                                                             &
137        ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, p, prho, pt, q, qc, ql, qr, qs, &
138               qsws, qswst, rho, sa, saswsb, saswst, shf, td_lsa_lpt, td_lsa_q,&
139               td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us, usws,    &
140               uswst, vsws, v, vg, vpt, vswst, w, w_subs, zw
141       
142    USE cloud_parameters,                                                      &
143        ONLY:   l_d_cp, prr, pt_d_t
144       
145    USE control_parameters,                                                    &
146        ONLY:   average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
147                dt_3d, g, humidity, icloud_scheme, kappa, large_scale_forcing, &
148                large_scale_subsidence, max_pr_user, message_string,           &
149                monotonic_adjustment, ocean,                                   &
150                passive_scalar, precipitation, simulated_time,                 &
151                use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
152                ws_scheme_mom, ws_scheme_sca
153       
154    USE cpulog,                                                                &
155        ONLY:   cpu_log, log_point
156       
157    USE grid_variables,                                                        &
158        ONLY:   ddx, ddy
159       
160    USE indices,                                                               &
161        ONLY:   ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums,      &
162                ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzb_diff_s_inner,        &
163                nzb_s_inner, nzt, nzt_diff
164       
165    USE kinds
166   
167    USE land_surface_model_mod,                                                &
168        ONLY:   dots_soil, ghf_eb, land_surface, m_soil, nzb_soil, nzt_soil,   &
169                qsws_eb, qsws_liq_eb, qsws_soil_eb, qsws_veg_eb, r_a, r_s,     &
170                shf_eb, t_soil
171
172    USE pegrid
173
174    USE radiation_model_mod,                                                   &
175        ONLY :  dots_rad, rad_net, rad_sw_in, radiation
176   
177    USE statistics
178
179    IMPLICIT NONE
180
181    INTEGER(iwp) ::  i                   !:
182    INTEGER(iwp) ::  j                   !:
183    INTEGER(iwp) ::  k                   !:
184    INTEGER(iwp) ::  nt                  !:
185    INTEGER(iwp) ::  omp_get_thread_num  !:
186    INTEGER(iwp) ::  sr                  !:
187    INTEGER(iwp) ::  tn                  !:
188   
189    LOGICAL ::  first  !:
190   
191    REAL(wp) ::  dptdz_threshold  !:
192    REAL(wp) ::  fac              !:
193    REAL(wp) ::  height           !:
194    REAL(wp) ::  pts              !:
195    REAL(wp) ::  sums_l_eper      !:
196    REAL(wp) ::  sums_l_etot      !:
197    REAL(wp) ::  ust              !:
198    REAL(wp) ::  ust2             !:
199    REAL(wp) ::  u2               !:
200    REAL(wp) ::  vst              !:
201    REAL(wp) ::  vst2             !:
202    REAL(wp) ::  v2               !:
203    REAL(wp) ::  w2               !:
204    REAL(wp) ::  z_i(2)           !:
205   
206    REAL(wp) ::  dptdz(nzb+1:nzt+1)    !:
207    REAL(wp) ::  sums_ll(nzb:nzt+1,2)  !:
208
209    CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
210
211    !$acc update host( km, kh, e, pt, qs, qsws, rif, shf, ts, u, usws, v, vsws, w )
212
213!
214!-- To be on the safe side, check whether flow_statistics has already been
215!-- called once after the current time step
216    IF ( flow_statistics_called )  THEN
217
218       message_string = 'flow_statistics is called two times within one ' // &
219                        'timestep'
220       CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
221
222    ENDIF
223
224!
225!-- Compute statistics for each (sub-)region
226    DO  sr = 0, statistic_regions
227
228!
229!--    Initialize (local) summation array
230       sums_l = 0.0_wp
231
232!
233!--    Store sums that have been computed in other subroutines in summation
234!--    array
235       sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
236!--    WARNING: next line still has to be adjusted for OpenMP
237       sums_l(:,21,0) = sums_wsts_bc_l(:,sr)  ! heat flux from advec_s_bc
238       sums_l(nzb+9,pr_palm,0)  = sums_divold_l(sr)  ! old divergence from pres
239       sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
240
241!
242!--    When calcuating horizontally-averaged total (resolved- plus subgrid-
243!--    scale) vertical fluxes and velocity variances by using commonly-
244!--    applied Reynolds-based methods ( e.g. <w'pt'> = (w-<w>)*(pt-<pt>) )
245!--    in combination with the 5th order advection scheme, pronounced
246!--    artificial kinks could be observed in the vertical profiles near the
247!--    surface. Please note: these kinks were not related to the model truth,
248!--    i.e. these kinks are just related to an evaluation problem.   
249!--    In order avoid these kinks, vertical fluxes and horizontal as well
250!--    vertical velocity variances are calculated directly within the advection
251!--    routines, according to the numerical discretization, to evaluate the
252!--    statistical quantities as they will appear within the prognostic
253!--    equations.
254!--    Copy the turbulent quantities, evaluated in the advection routines to
255!--    the local array sums_l() for further computations.
256       IF ( ws_scheme_mom .AND. sr == 0 )  THEN
257
258!
259!--       According to the Neumann bc for the horizontal velocity components,
260!--       the corresponding fluxes has to satisfiy the same bc.
261          IF ( ocean )  THEN
262             sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
263             sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)
264          ENDIF
265
266          DO  i = 0, threads_per_task-1
267!
268!--          Swap the turbulent quantities evaluated in advec_ws.
269             sums_l(:,13,i) = sums_wsus_ws_l(:,i)       ! w*u*
270             sums_l(:,15,i) = sums_wsvs_ws_l(:,i)       ! w*v*
271             sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2
272             sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2
273             sums_l(:,32,i) = sums_ws2_ws_l(:,i)        ! w*2
274             sums_l(:,34,i) = sums_l(:,34,i) + 0.5_wp *                        & 
275                              ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +      &
276                                sums_ws2_ws_l(:,i) )    ! e*
277             DO  k = nzb, nzt
278                sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5_wp * ( &
279                                                      sums_us2_ws_l(k,i) +     &
280                                                      sums_vs2_ws_l(k,i) +     &
281                                                      sums_ws2_ws_l(k,i) )
282             ENDDO
283          ENDDO
284
285       ENDIF
286
287       IF ( ws_scheme_sca .AND. .NOT. monotonic_adjustment                    &
288           .AND. sr == 0 )  THEN
289
290          DO  i = 0, threads_per_task-1
291             sums_l(:,17,i) = sums_wspts_ws_l(:,i)      ! w*pt* from advec_s_ws
292             IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
293             IF ( humidity .OR. passive_scalar ) sums_l(:,49,i) =              &
294                                                   sums_wsqs_ws_l(:,i) !w*q*
295          ENDDO
296
297       ENDIF
298!
299!--    Horizontally averaged profiles of horizontal velocities and temperature.
300!--    They must have been computed before, because they are already required
301!--    for other horizontal averages.
302       tn = 0
303
304       !$OMP PARALLEL PRIVATE( i, j, k, tn )
305#if defined( __intel_openmp_bug )
306       tn = omp_get_thread_num()
307#else
308!$     tn = omp_get_thread_num()
309#endif
310
311       !$OMP DO
312       DO  i = nxl, nxr
313          DO  j =  nys, nyn
314             DO  k = nzb_s_inner(j,i), nzt+1
315                sums_l(k,1,tn)  = sums_l(k,1,tn)  + u(k,j,i)  * rmask(j,i,sr)
316                sums_l(k,2,tn)  = sums_l(k,2,tn)  + v(k,j,i)  * rmask(j,i,sr)
317                sums_l(k,4,tn)  = sums_l(k,4,tn)  + pt(k,j,i) * rmask(j,i,sr)
318             ENDDO
319          ENDDO
320       ENDDO
321
322!
323!--    Horizontally averaged profile of salinity
324       IF ( ocean )  THEN
325          !$OMP DO
326          DO  i = nxl, nxr
327             DO  j =  nys, nyn
328                DO  k = nzb_s_inner(j,i), nzt+1
329                   sums_l(k,23,tn)  = sums_l(k,23,tn) + &
330                                      sa(k,j,i) * rmask(j,i,sr)
331                ENDDO
332             ENDDO
333          ENDDO
334       ENDIF
335
336!
337!--    Horizontally averaged profiles of virtual potential temperature,
338!--    total water content, specific humidity and liquid water potential
339!--    temperature
340       IF ( humidity )  THEN
341          !$OMP DO
342          DO  i = nxl, nxr
343             DO  j =  nys, nyn
344                DO  k = nzb_s_inner(j,i), nzt+1
345                   sums_l(k,44,tn)  = sums_l(k,44,tn) + &
346                                      vpt(k,j,i) * rmask(j,i,sr)
347                   sums_l(k,41,tn)  = sums_l(k,41,tn) + &
348                                      q(k,j,i) * rmask(j,i,sr)
349                ENDDO
350             ENDDO
351          ENDDO
352          IF ( cloud_physics )  THEN
353             !$OMP DO
354             DO  i = nxl, nxr
355                DO  j =  nys, nyn
356                   DO  k = nzb_s_inner(j,i), nzt+1
357                      sums_l(k,42,tn) = sums_l(k,42,tn) + &
358                                      ( q(k,j,i) - ql(k,j,i) ) * rmask(j,i,sr)
359                      sums_l(k,43,tn) = sums_l(k,43,tn) + ( &
360                                      pt(k,j,i) + l_d_cp*pt_d_t(k) * ql(k,j,i) &
361                                                          ) * rmask(j,i,sr)
362                   ENDDO
363                ENDDO
364             ENDDO
365          ENDIF
366       ENDIF
367
368!
369!--    Horizontally averaged profiles of passive scalar
370       IF ( passive_scalar )  THEN
371          !$OMP DO
372          DO  i = nxl, nxr
373             DO  j =  nys, nyn
374                DO  k = nzb_s_inner(j,i), nzt+1
375                   sums_l(k,41,tn)  = sums_l(k,41,tn) + q(k,j,i) * rmask(j,i,sr)
376                ENDDO
377             ENDDO
378          ENDDO
379       ENDIF
380       !$OMP END PARALLEL
381!
382!--    Summation of thread sums
383       IF ( threads_per_task > 1 )  THEN
384          DO  i = 1, threads_per_task-1
385             sums_l(:,1,0) = sums_l(:,1,0) + sums_l(:,1,i)
386             sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i)
387             sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i)
388             IF ( ocean )  THEN
389                sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)
390             ENDIF
391             IF ( humidity )  THEN
392                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
393                sums_l(:,44,0) = sums_l(:,44,0) + sums_l(:,44,i)
394                IF ( cloud_physics )  THEN
395                   sums_l(:,42,0) = sums_l(:,42,0) + sums_l(:,42,i)
396                   sums_l(:,43,0) = sums_l(:,43,0) + sums_l(:,43,i)
397                ENDIF
398             ENDIF
399             IF ( passive_scalar )  THEN
400                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
401             ENDIF
402          ENDDO
403       ENDIF
404
405#if defined( __parallel )
406!
407!--    Compute total sum from local sums
408       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
409       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL,  &
410                           MPI_SUM, comm2d, ierr )
411       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
412       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL,  &
413                           MPI_SUM, comm2d, ierr )
414       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
415       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL,  &
416                           MPI_SUM, comm2d, ierr )
417       IF ( ocean )  THEN
418          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
419          CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb,       &
420                              MPI_REAL, MPI_SUM, comm2d, ierr )
421       ENDIF
422       IF ( humidity ) THEN
423          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
424          CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb,       &
425                              MPI_REAL, MPI_SUM, comm2d, ierr )
426          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
427          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb,       &
428                              MPI_REAL, MPI_SUM, comm2d, ierr )
429          IF ( cloud_physics ) THEN
430             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
431             CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb,    &
432                                 MPI_REAL, MPI_SUM, comm2d, ierr )
433             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
434             CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb,    &
435                                 MPI_REAL, MPI_SUM, comm2d, ierr )
436          ENDIF
437       ENDIF
438
439       IF ( passive_scalar )  THEN
440          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
441          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb,       &
442                              MPI_REAL, MPI_SUM, comm2d, ierr )
443       ENDIF
444#else
445       sums(:,1) = sums_l(:,1,0)
446       sums(:,2) = sums_l(:,2,0)
447       sums(:,4) = sums_l(:,4,0)
448       IF ( ocean )  sums(:,23) = sums_l(:,23,0)
449       IF ( humidity ) THEN
450          sums(:,44) = sums_l(:,44,0)
451          sums(:,41) = sums_l(:,41,0)
452          IF ( cloud_physics ) THEN
453             sums(:,42) = sums_l(:,42,0)
454             sums(:,43) = sums_l(:,43,0)
455          ENDIF
456       ENDIF
457       IF ( passive_scalar )  sums(:,41) = sums_l(:,41,0)
458#endif
459
460!
461!--    Final values are obtained by division by the total number of grid points
462!--    used for summation. After that store profiles.
463       sums(:,1) = sums(:,1) / ngp_2dh(sr)
464       sums(:,2) = sums(:,2) / ngp_2dh(sr)
465       sums(:,4) = sums(:,4) / ngp_2dh_s_inner(:,sr)
466       hom(:,1,1,sr) = sums(:,1)             ! u
467       hom(:,1,2,sr) = sums(:,2)             ! v
468       hom(:,1,4,sr) = sums(:,4)             ! pt
469
470
471!
472!--    Salinity
473       IF ( ocean )  THEN
474          sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr)
475          hom(:,1,23,sr) = sums(:,23)             ! sa
476       ENDIF
477
478!
479!--    Humidity and cloud parameters
480       IF ( humidity ) THEN
481          sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr)
482          sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
483          hom(:,1,44,sr) = sums(:,44)             ! vpt
484          hom(:,1,41,sr) = sums(:,41)             ! qv (q)
485          IF ( cloud_physics ) THEN
486             sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr)
487             sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr)
488             hom(:,1,42,sr) = sums(:,42)             ! qv
489             hom(:,1,43,sr) = sums(:,43)             ! pt
490          ENDIF
491       ENDIF
492
493!
494!--    Passive scalar
495       IF ( passive_scalar )  hom(:,1,41,sr) = sums(:,41) /                    &
496            ngp_2dh_s_inner(:,sr)                    ! s (q)
497
498!
499!--    Horizontally averaged profiles of the remaining prognostic variables,
500!--    variances, the total and the perturbation energy (single values in last
501!--    column of sums_l) and some diagnostic quantities.
502!--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
503!--    ----  speaking the following k-loop would have to be split up and
504!--          rearranged according to the staggered grid.
505!--          However, this implies no error since staggered velocity components
506!--          are zero at the walls and inside buildings.
507       tn = 0
508#if defined( __intel_openmp_bug )
509       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, &
510       !$OMP                    tn, ust, ust2, u2, vst, vst2, v2, w2 )
511       tn = omp_get_thread_num()
512#else
513       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2, w2 )
514!$     tn = omp_get_thread_num()
515#endif
516       !$OMP DO
517       DO  i = nxl, nxr
518          DO  j =  nys, nyn
519             sums_l_etot = 0.0_wp
520             DO  k = nzb_s_inner(j,i), nzt+1
521!
522!--             Prognostic and diagnostic variables
523                sums_l(k,3,tn)  = sums_l(k,3,tn)  + w(k,j,i)  * rmask(j,i,sr)
524                sums_l(k,8,tn)  = sums_l(k,8,tn)  + e(k,j,i)  * rmask(j,i,sr)
525                sums_l(k,9,tn)  = sums_l(k,9,tn)  + km(k,j,i) * rmask(j,i,sr)
526                sums_l(k,10,tn) = sums_l(k,10,tn) + kh(k,j,i) * rmask(j,i,sr)
527                sums_l(k,40,tn) = sums_l(k,40,tn) + p(k,j,i)
528
529                sums_l(k,33,tn) = sums_l(k,33,tn) + &
530                                  ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr)
531
532                IF ( humidity )  THEN
533                   sums_l(k,70,tn) = sums_l(k,70,tn) + &
534                                  ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr)
535                ENDIF
536
537!
538!--             Higher moments
539!--             (Computation of the skewness of w further below)
540                sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i)**3 * rmask(j,i,sr)
541
542                sums_l_etot  = sums_l_etot + &
543                                        0.5_wp * ( u(k,j,i)**2 + v(k,j,i)**2 + &
544                                        w(k,j,i)**2 ) * rmask(j,i,sr)
545             ENDDO
546!
547!--          Total and perturbation energy for the total domain (being
548!--          collected in the last column of sums_l). Summation of these
549!--          quantities is seperated from the previous loop in order to
550!--          allow vectorization of that loop.
551             sums_l(nzb+4,pr_palm,tn) = sums_l(nzb+4,pr_palm,tn) + sums_l_etot
552!
553!--          2D-arrays (being collected in the last column of sums_l)
554             sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +               &
555                                        us(j,i)   * rmask(j,i,sr)
556             sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +             &
557                                        usws(j,i) * rmask(j,i,sr)
558             sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +             &
559                                        vsws(j,i) * rmask(j,i,sr)
560             sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +             &
561                                        ts(j,i)   * rmask(j,i,sr)
562             IF ( humidity )  THEN
563                sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +        &
564                                            qs(j,i)   * rmask(j,i,sr)
565             ENDIF
566          ENDDO
567       ENDDO
568
569!
570!--    Computation of statistics when ws-scheme is not used. Else these
571!--    quantities are evaluated in the advection routines.
572       IF ( .NOT. ws_scheme_mom .OR. sr /= 0 )  THEN
573          !$OMP DO
574          DO  i = nxl, nxr
575             DO  j =  nys, nyn
576                sums_l_eper = 0.0_wp
577                DO  k = nzb_s_inner(j,i), nzt+1
578                   u2   = u(k,j,i)**2
579                   v2   = v(k,j,i)**2
580                   w2   = w(k,j,i)**2
581                   ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
582                   vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
583
584                   sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr)
585                   sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr)
586                   sums_l(k,32,tn) = sums_l(k,32,tn) + w2   * rmask(j,i,sr)
587!
588!--             Perturbation energy
589
590                   sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5_wp *                &
591                                  ( ust2 + vst2 + w2 ) * rmask(j,i,sr)
592                   sums_l_eper     = sums_l_eper +                             &
593                                     0.5_wp * ( ust2+vst2+w2 ) * rmask(j,i,sr)
594
595                ENDDO
596                sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn)            &
597                     + sums_l_eper
598             ENDDO
599          ENDDO
600       ENDIF
601
602!
603!--    Horizontally averaged profiles of the vertical fluxes
604
605       !$OMP DO
606       DO  i = nxl, nxr
607          DO  j = nys, nyn
608!
609!--          Subgridscale fluxes (without Prandtl layer from k=nzb,
610!--          oterwise from k=nzb+1)
611!--          NOTE: for simplicity, nzb_diff_s_inner is used below, although
612!--          ----  strictly speaking the following k-loop would have to be
613!--                split up according to the staggered grid.
614!--                However, this implies no error since staggered velocity
615!--                components are zero at the walls and inside buildings.
616
617             DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
618!
619!--             Momentum flux w"u"
620                sums_l(k,12,tn) = sums_l(k,12,tn) - 0.25_wp * (                &
621                               km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
622                                                           ) * (               &
623                                   ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
624                                 + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
625                                                               ) * rmask(j,i,sr)
626!
627!--             Momentum flux w"v"
628                sums_l(k,14,tn) = sums_l(k,14,tn) - 0.25_wp * (                &
629                               km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
630                                                           ) * (               &
631                                   ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
632                                 + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
633                                                               ) * rmask(j,i,sr)
634!
635!--             Heat flux w"pt"
636                sums_l(k,16,tn) = sums_l(k,16,tn)                              &
637                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
638                                               * ( pt(k+1,j,i) - pt(k,j,i) )   &
639                                               * ddzu(k+1) * rmask(j,i,sr)
640
641
642!
643!--             Salinity flux w"sa"
644                IF ( ocean )  THEN
645                   sums_l(k,65,tn) = sums_l(k,65,tn)                           &
646                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
647                                               * ( sa(k+1,j,i) - sa(k,j,i) )   &
648                                               * ddzu(k+1) * rmask(j,i,sr)
649                ENDIF
650
651!
652!--             Buoyancy flux, water flux (humidity flux) w"q"
653                IF ( humidity ) THEN
654                   sums_l(k,45,tn) = sums_l(k,45,tn)                           &
655                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
656                                               * ( vpt(k+1,j,i) - vpt(k,j,i) ) &
657                                               * ddzu(k+1) * rmask(j,i,sr)
658                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
659                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
660                                               * ( q(k+1,j,i) - q(k,j,i) )     &
661                                               * ddzu(k+1) * rmask(j,i,sr)
662
663                   IF ( cloud_physics ) THEN
664                      sums_l(k,51,tn) = sums_l(k,51,tn)                        &
665                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
666                                               * ( ( q(k+1,j,i) - ql(k+1,j,i) )&
667                                                - ( q(k,j,i) - ql(k,j,i) ) )   &
668                                               * ddzu(k+1) * rmask(j,i,sr) 
669                   ENDIF
670                ENDIF
671
672!
673!--             Passive scalar flux
674                IF ( passive_scalar )  THEN
675                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
676                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
677                                               * ( q(k+1,j,i) - q(k,j,i) )     &
678                                               * ddzu(k+1) * rmask(j,i,sr)
679                ENDIF
680
681             ENDDO
682
683!
684!--          Subgridscale fluxes in the Prandtl layer
685             IF ( use_surface_fluxes )  THEN
686                sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
687                                    usws(j,i) * rmask(j,i,sr)     ! w"u"
688                sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
689                                    vsws(j,i) * rmask(j,i,sr)     ! w"v"
690                sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
691                                    shf(j,i)  * rmask(j,i,sr)     ! w"pt"
692                sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
693                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
694                sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
695                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
696                IF ( ocean )  THEN
697                   sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
698                                       saswsb(j,i) * rmask(j,i,sr)  ! w"sa"
699                ENDIF
700                IF ( humidity )  THEN
701                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
702                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
703                   sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                   &
704                                       ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) *     &
705                                       shf(j,i) + 0.61_wp * pt(nzb,j,i) *      &
706                                                  qsws(j,i) )
707                   IF ( cloud_droplets )  THEN
708                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                &
709                                         ( 1.0_wp + 0.61_wp * q(nzb,j,i) -     &
710                                           ql(nzb,j,i) ) * shf(j,i) +          &
711                                           0.61_wp * pt(nzb,j,i) * qsws(j,i) )
712                   ENDIF
713                   IF ( cloud_physics )  THEN
714!
715!--                   Formula does not work if ql(nzb) /= 0.0
716                      sums_l(nzb,51,tn) = sums_l(nzb,51,tn) + &   ! w"q" (w"qv")
717                                          qsws(j,i) * rmask(j,i,sr)
718                   ENDIF
719                ENDIF
720                IF ( passive_scalar )  THEN
721                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
722                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
723                ENDIF
724             ENDIF
725
726             IF ( land_surface )  THEN
727                sums_l(nzb,93,tn)  = sums_l(nzb,93,tn) + ghf_eb(j,i)
728                sums_l(nzb,94,tn)  = sums_l(nzb,94,tn) + shf_eb(j,i)
729                sums_l(nzb,95,tn)  = sums_l(nzb,95,tn) + qsws_eb(j,i)
730                sums_l(nzb,96,tn)  = sums_l(nzb,96,tn) + qsws_liq_eb(j,i)
731                sums_l(nzb,97,tn)  = sums_l(nzb,97,tn) + qsws_soil_eb(j,i)
732                sums_l(nzb,98,tn)  = sums_l(nzb,98,tn) + qsws_veg_eb(j,i)
733                sums_l(nzb,99,tn)  = sums_l(nzb,99,tn) + r_a(j,i)
734                sums_l(nzb,100,tn) = sums_l(nzb,100,tn)+ r_s(j,i)
735             ENDIF
736
737             IF ( radiation )  THEN
738                sums_l(nzb,101,tn)  = sums_l(nzb,101,tn)  + rad_net(j,i)
739                sums_l(nzb,102,tn) = sums_l(nzb,102,tn) + rad_sw_in(j,i)
740             ENDIF
741
742!
743!--          Subgridscale fluxes at the top surface
744             IF ( use_top_fluxes )  THEN
745                sums_l(nzt:nzt+1,12,tn) = sums_l(nzt:nzt+1,12,tn) + &
746                                    uswst(j,i) * rmask(j,i,sr)    ! w"u"
747                sums_l(nzt:nzt+1,14,tn) = sums_l(nzt:nzt+1,14,tn) + &
748                                    vswst(j,i) * rmask(j,i,sr)    ! w"v"
749                sums_l(nzt:nzt+1,16,tn) = sums_l(nzt:nzt+1,16,tn) + &
750                                    tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
751                sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + &
752                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
753                sums_l(nzt:nzt+1,61,tn) = sums_l(nzt:nzt+1,61,tn) + &
754                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
755
756                IF ( ocean )  THEN
757                   sums_l(nzt,65,tn) = sums_l(nzt,65,tn) + &
758                                       saswst(j,i) * rmask(j,i,sr)  ! w"sa"
759                ENDIF
760                IF ( humidity )  THEN
761                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) +                     &
762                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
763                   sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                   &
764                                       ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) *     &
765                                       tswst(j,i) + 0.61_wp * pt(nzt,j,i) *    &
766                                                             qswst(j,i) )
767                   IF ( cloud_droplets )  THEN
768                      sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                &
769                                          ( 1.0_wp + 0.61_wp * q(nzt,j,i) -    &
770                                            ql(nzt,j,i) ) * tswst(j,i) +       &
771                                            0.61_wp * pt(nzt,j,i) * qswst(j,i) )
772                   ENDIF
773                   IF ( cloud_physics )  THEN
774!
775!--                   Formula does not work if ql(nzb) /= 0.0
776                      sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + &   ! w"q" (w"qv")
777                                          qswst(j,i) * rmask(j,i,sr)
778                   ENDIF
779                ENDIF
780                IF ( passive_scalar )  THEN
781                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
782                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
783                ENDIF
784             ENDIF
785
786!
787!--          Resolved fluxes (can be computed for all horizontal points)
788!--          NOTE: for simplicity, nzb_s_inner is used below, although strictly
789!--          ----  speaking the following k-loop would have to be split up and
790!--                rearranged according to the staggered grid.
791             DO  k = nzb_s_inner(j,i), nzt
792                ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +                  &
793                                 u(k+1,j,i) - hom(k+1,1,1,sr) )
794                vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +                  &
795                                 v(k+1,j,i) - hom(k+1,1,2,sr) )
796                pts = 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) +                 &
797                                 pt(k+1,j,i) - hom(k+1,1,4,sr) )
798
799!--             Higher moments
800                sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 *        &
801                                                    rmask(j,i,sr)
802                sums_l(k,36,tn) = sums_l(k,36,tn) + pts**2 * w(k,j,i) *        &
803                                                    rmask(j,i,sr)
804
805!
806!--             Salinity flux and density (density does not belong to here,
807!--             but so far there is no other suitable place to calculate)
808                IF ( ocean )  THEN
809                   IF( .NOT. ws_scheme_sca .OR. monotonic_adjustment .OR.      &
810                             sr /= 0 )  THEN
811                      pts = 0.5_wp * ( sa(k,j,i)   - hom(k,1,23,sr) +          &
812                                       sa(k+1,j,i) - hom(k+1,1,23,sr) )
813                      sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) *     &
814                                        rmask(j,i,sr)
815                   ENDIF
816                   sums_l(k,64,tn) = sums_l(k,64,tn) + rho(k,j,i) *            &
817                                                       rmask(j,i,sr)
818                   sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) *           &
819                                                       rmask(j,i,sr)
820                ENDIF
821
822!
823!--             Buoyancy flux, water flux, humidity flux, liquid water
824!--             content, rain drop concentration and rain water content
825                IF ( humidity )  THEN
826                   IF ( cloud_physics .OR. cloud_droplets )  THEN
827                      pts = 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +         &
828                                    vpt(k+1,j,i) - hom(k+1,1,44,sr) )
829                      sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *     &
830                                                          rmask(j,i,sr)
831                      IF ( .NOT. cloud_droplets )  THEN
832                         pts = 0.5_wp *                                        &
833                              ( ( q(k,j,i) - ql(k,j,i) ) -                     &
834                              hom(k,1,42,sr) +                                 &
835                              ( q(k+1,j,i) - ql(k+1,j,i) ) -                   &
836                              hom(k+1,1,42,sr) )
837                         sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) *  &
838                                                             rmask(j,i,sr)
839                         IF ( icloud_scheme == 0  )  THEN
840                            sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) *    &
841                                                                rmask(j,i,sr)
842                            sums_l(k,75,tn) = sums_l(k,75,tn) + qc(k,j,i) *    &
843                                                                rmask(j,i,sr)
844                            IF ( precipitation )  THEN
845                               sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) * &
846                                                                   rmask(j,i,sr)
847                               sums_l(k,74,tn) = sums_l(k,74,tn) + qr(k,j,i) * &
848                                                                   rmask(j,i,sr)
849                               sums_l(k,76,tn) = sums_l(k,76,tn) + prr(k,j,i) *&
850                                                                   rmask(j,i,sr)
851                            ENDIF
852                         ELSE
853                            sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) *    &
854                                                                rmask(j,i,sr)
855                         ENDIF
856                      ELSE
857                         sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) *       &
858                                                             rmask(j,i,sr)
859                      ENDIF
860                   ELSE
861                      IF( .NOT. ws_scheme_sca .OR. monotonic_adjustment .OR.   &
862                                sr /= 0 )  THEN
863                         pts = 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +      &
864                                          vpt(k+1,j,i) - hom(k+1,1,44,sr) )
865                         sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *  &
866                                                             rmask(j,i,sr)
867                      ELSE IF ( ws_scheme_sca .AND. .NOT. monotonic_adjustment &
868                                .AND. sr == 0 )  THEN
869                         sums_l(k,46,tn) = ( 1.0_wp + 0.61_wp *                & 
870                                             hom(k,1,41,sr) ) *                &
871                                           sums_l(k,17,tn) +                   &
872                                           0.61_wp * hom(k,1,4,sr) *           &
873                                           sums_l(k,49,tn)
874                      END IF
875                   END IF
876                ENDIF
877!
878!--             Passive scalar flux
879                IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca                &
880                     .OR. monotonic_adjustment .OR. sr /= 0 ) )  THEN
881                   pts = 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +              &
882                                    q(k+1,j,i) - hom(k+1,1,41,sr) )
883                   sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *        &
884                                                       rmask(j,i,sr)
885                ENDIF
886
887!
888!--             Energy flux w*e*
889!--             has to be adjusted
890                sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5_wp *        &
891                                             ( ust**2 + vst**2 + w(k,j,i)**2 ) &
892                                             * rmask(j,i,sr)
893             ENDDO
894          ENDDO
895       ENDDO
896!
897!--    For speed optimization fluxes which have been computed in part directly
898!--    inside the WS advection routines are treated seperatly
899!--    Momentum fluxes first:
900       IF ( .NOT. ws_scheme_mom .OR. sr /= 0  )  THEN
901         !$OMP DO
902         DO  i = nxl, nxr
903            DO  j = nys, nyn
904               DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
905                  ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +                &
906                                   u(k+1,j,i) - hom(k+1,1,1,sr) )
907                  vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +                &
908                                   v(k+1,j,i) - hom(k+1,1,2,sr) )
909!
910!--               Momentum flux w*u*
911                  sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5_wp *                 &
912                                                    ( w(k,j,i-1) + w(k,j,i) )  &
913                                                    * ust * rmask(j,i,sr)
914!
915!--               Momentum flux w*v*
916                  sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5_wp *                 &
917                                                    ( w(k,j-1,i) + w(k,j,i) )  &
918                                                    * vst * rmask(j,i,sr)
919               ENDDO
920            ENDDO
921         ENDDO
922
923       ENDIF
924       IF ( .NOT. ws_scheme_sca .OR. monotonic_adjustment .OR. sr /= 0 )  THEN
925         !$OMP DO
926         DO  i = nxl, nxr
927            DO  j = nys, nyn
928               DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
929!
930!--               Vertical heat flux
931                  sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5_wp *                 &
932                           ( pt(k,j,i)   - hom(k,1,4,sr) +                     &
933                             pt(k+1,j,i) - hom(k+1,1,4,sr) )                   &
934                           * w(k,j,i) * rmask(j,i,sr)
935                  IF ( humidity )  THEN
936                     pts = 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +            &
937                                      q(k+1,j,i) - hom(k+1,1,41,sr) )
938                     sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *      &
939                                       rmask(j,i,sr)
940                  ENDIF
941               ENDDO
942            ENDDO
943         ENDDO
944
945       ENDIF
946
947!
948!--    Density at top follows Neumann condition
949       IF ( ocean )  THEN
950          sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
951          sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
952       ENDIF
953
954!
955!--    Divergence of vertical flux of resolved scale energy and pressure
956!--    fluctuations as well as flux of pressure fluctuation itself (68).
957!--    First calculate the products, then the divergence.
958!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
959       IF ( hom(nzb+1,2,55,0) /= 0.0_wp  .OR.  hom(nzb+1,2,68,0) /= 0.0_wp )  THEN
960
961          sums_ll = 0.0_wp  ! local array
962
963          !$OMP DO
964          DO  i = nxl, nxr
965             DO  j = nys, nyn
966                DO  k = nzb_s_inner(j,i)+1, nzt
967
968                   sums_ll(k,1) = sums_ll(k,1) + 0.5_wp * w(k,j,i) * (         &
969                  ( 0.25_wp * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1)    &
970                            - 0.5_wp * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) )     &
971                              ) )**2                                           &
972                + ( 0.25_wp * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i)    &
973                            - 0.5_wp * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) )     &
974                              ) )**2                                           &
975                + w(k,j,i)**2                                        )
976
977                   sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i)             &
978                                               * ( p(k,j,i) + p(k+1,j,i) )
979
980                ENDDO
981             ENDDO
982          ENDDO
983          sums_ll(0,1)     = 0.0_wp    ! because w is zero at the bottom
984          sums_ll(nzt+1,1) = 0.0_wp
985          sums_ll(0,2)     = 0.0_wp
986          sums_ll(nzt+1,2) = 0.0_wp
987
988          DO  k = nzb+1, nzt
989             sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
990             sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
991             sums_l(k,68,tn) = sums_ll(k,2)
992          ENDDO
993          sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
994          sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
995          sums_l(nzb,68,tn) = 0.0_wp    ! because w* = 0 at nzb
996
997       ENDIF
998
999!
1000!--    Divergence of vertical flux of SGS TKE and the flux itself (69)
1001       IF ( hom(nzb+1,2,57,0) /= 0.0_wp  .OR.  hom(nzb+1,2,69,0) /= 0.0_wp )  THEN
1002
1003          !$OMP DO
1004          DO  i = nxl, nxr
1005             DO  j = nys, nyn
1006                DO  k = nzb_s_inner(j,i)+1, nzt
1007
1008                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5_wp * (              &
1009                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
1010                 - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
1011                                                                ) * ddzw(k)
1012
1013                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5_wp * (              &
1014                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
1015                                                                )
1016
1017                ENDDO
1018             ENDDO
1019          ENDDO
1020          sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
1021          sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
1022
1023       ENDIF
1024
1025!
1026!--    Horizontal heat fluxes (subgrid, resolved, total).
1027!--    Do it only, if profiles shall be plotted.
1028       IF ( hom(nzb+1,2,58,0) /= 0.0_wp ) THEN
1029
1030          !$OMP DO
1031          DO  i = nxl, nxr
1032             DO  j = nys, nyn
1033                DO  k = nzb_s_inner(j,i)+1, nzt
1034!
1035!--                Subgrid horizontal heat fluxes u"pt", v"pt"
1036                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5_wp *                &
1037                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
1038                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
1039                                                 * ddx * rmask(j,i,sr)
1040                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
1041                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
1042                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
1043                                                 * ddy * rmask(j,i,sr)
1044!
1045!--                Resolved horizontal heat fluxes u*pt*, v*pt*
1046                   sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
1047                                                  ( u(k,j,i) - hom(k,1,1,sr) ) &
1048                                    * 0.5_wp * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
1049                                                 pt(k,j,i)   - hom(k,1,4,sr) )
1050                   pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
1051                                    pt(k,j,i)   - hom(k,1,4,sr) )
1052                   sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
1053                                                  ( v(k,j,i) - hom(k,1,2,sr) ) &
1054                                    * 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
1055                                                 pt(k,j,i)   - hom(k,1,4,sr) )
1056                ENDDO
1057             ENDDO
1058          ENDDO
1059!
1060!--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
1061          sums_l(nzb,58,tn) = 0.0_wp
1062          sums_l(nzb,59,tn) = 0.0_wp
1063          sums_l(nzb,60,tn) = 0.0_wp
1064          sums_l(nzb,61,tn) = 0.0_wp
1065          sums_l(nzb,62,tn) = 0.0_wp
1066          sums_l(nzb,63,tn) = 0.0_wp
1067
1068       ENDIF
1069
1070!
1071!--    Collect current large scale advection and subsidence tendencies for
1072!--    data output
1073       IF ( large_scale_forcing  .AND.  ( simulated_time .GT. 0.0_wp ) )  THEN
1074!
1075!--       Interpolation in time of LSF_DATA
1076          nt = 1
1077          DO WHILE ( simulated_time - dt_3d > time_vert(nt) )
1078             nt = nt + 1
1079          ENDDO
1080          IF ( simulated_time - dt_3d /= time_vert(nt) )  THEN
1081            nt = nt - 1
1082          ENDIF
1083
1084          fac = ( simulated_time - dt_3d - time_vert(nt) )                     &
1085                / ( time_vert(nt+1)-time_vert(nt) )
1086
1087
1088          DO  k = nzb, nzt
1089             sums_ls_l(k,0) = td_lsa_lpt(k,nt)                                 &
1090                              + fac * ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )
1091             sums_ls_l(k,1) = td_lsa_q(k,nt)                                   &
1092                              + fac * ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) )
1093          ENDDO
1094
1095          sums_ls_l(nzt+1,0) = sums_ls_l(nzt,0)
1096          sums_ls_l(nzt+1,1) = sums_ls_l(nzt,1)
1097
1098          IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
1099
1100             DO  k = nzb, nzt
1101                sums_ls_l(k,2) = td_sub_lpt(k,nt) + fac *                      &
1102                                 ( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )
1103                sums_ls_l(k,3) = td_sub_q(k,nt) + fac *                        &
1104                                 ( td_sub_q(k,nt+1) - td_sub_q(k,nt) )
1105             ENDDO
1106
1107             sums_ls_l(nzt+1,2) = sums_ls_l(nzt,2)
1108             sums_ls_l(nzt+1,3) = sums_ls_l(nzt,3)
1109
1110          ENDIF
1111
1112       ENDIF
1113
1114
1115       IF ( land_surface )  THEN
1116          !$OMP DO
1117          DO  i = nxl, nxr
1118             DO  j =  nys, nyn
1119                DO  k = nzb_soil, nzt_soil
1120                   sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil(k,j,i) * rmask(j,i,sr)
1121                   sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil(k,j,i) * rmask(j,i,sr)
1122                ENDDO
1123             ENDDO
1124          ENDDO
1125       ENDIF
1126       
1127
1128!
1129!--    Calculate the user-defined profiles
1130       CALL user_statistics( 'profiles', sr, tn )
1131       !$OMP END PARALLEL
1132
1133!
1134!--    Summation of thread sums
1135       IF ( threads_per_task > 1 )  THEN
1136          DO  i = 1, threads_per_task-1
1137             sums_l(:,3,0)          = sums_l(:,3,0) + sums_l(:,3,i)
1138             sums_l(:,4:40,0)       = sums_l(:,4:40,0) + sums_l(:,4:40,i)
1139             sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
1140                                      sums_l(:,45:pr_palm,i)
1141             IF ( max_pr_user > 0 )  THEN
1142                sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
1143                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
1144                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
1145             ENDIF
1146          ENDDO
1147       ENDIF
1148
1149#if defined( __parallel )
1150
1151!
1152!--    Compute total sum from local sums
1153       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1154       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL,   &
1155                           MPI_SUM, comm2d, ierr )
1156       IF ( large_scale_forcing )  THEN
1157          CALL MPI_ALLREDUCE( sums_ls_l(nzb,2), sums(nzb,83), ngp_sums_ls,     &
1158                              MPI_REAL, MPI_SUM, comm2d, ierr )
1159       ENDIF
1160#else
1161       sums = sums_l(:,:,0)
1162       IF ( large_scale_forcing )  THEN
1163          sums(:,81:88) = sums_ls_l
1164       ENDIF
1165#endif
1166
1167!
1168!--    Final values are obtained by division by the total number of grid points
1169!--    used for summation. After that store profiles.
1170!--    Profiles:
1171       DO  k = nzb, nzt+1
1172          sums(k,3)               = sums(k,3)           / ngp_2dh(sr)
1173          sums(k,8:11)            = sums(k,8:11)        / ngp_2dh_s_inner(k,sr)
1174          sums(k,12:22)           = sums(k,12:22)       / ngp_2dh(sr)
1175          sums(k,23:29)           = sums(k,23:29)       / ngp_2dh_s_inner(k,sr)
1176          sums(k,30:32)           = sums(k,30:32)       / ngp_2dh(sr)
1177          sums(k,33:34)           = sums(k,33:34)       / ngp_2dh_s_inner(k,sr)
1178          sums(k,35:39)           = sums(k,35:39)       / ngp_2dh(sr)
1179          sums(k,40)              = sums(k,40)          / ngp_2dh_s_inner(k,sr)
1180          sums(k,45:53)           = sums(k,45:53)       / ngp_2dh(sr)
1181          sums(k,54)              = sums(k,54)          / ngp_2dh_s_inner(k,sr)
1182          sums(k,55:63)           = sums(k,55:63)       / ngp_2dh(sr)
1183          sums(k,64)              = sums(k,64)          / ngp_2dh_s_inner(k,sr)
1184          sums(k,65:69)           = sums(k,65:69)       / ngp_2dh(sr)
1185          sums(k,70:80)           = sums(k,70:80)       / ngp_2dh_s_inner(k,sr)
1186          sums(k,81:88)           = sums(k,81:88)       / ngp_2dh(sr)
1187          sums(k,89:100)           = sums(k,89:100)     / ngp_2dh(sr)
1188          sums(k,101:pr_palm-2)    = sums(k,101:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
1189       ENDDO
1190
1191!--    Upstream-parts
1192       sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
1193!--    u* and so on
1194!--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
1195!--    size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
1196!--    above the topography, they are being divided by ngp_2dh(sr)
1197       sums(nzb:nzb+3,pr_palm)    = sums(nzb:nzb+3,pr_palm)    / &
1198                                    ngp_2dh(sr)
1199       sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)       / &    ! qs
1200                                    ngp_2dh(sr)
1201!--    eges, e*
1202       sums(nzb+4:nzb+5,pr_palm)  = sums(nzb+4:nzb+5,pr_palm)  / &
1203                                    ngp_3d(sr)
1204!--    Old and new divergence
1205       sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
1206                                    ngp_3d_inner(sr)
1207
1208!--    User-defined profiles
1209       IF ( max_pr_user > 0 )  THEN
1210          DO  k = nzb, nzt+1
1211             sums(k,pr_palm+1:pr_palm+max_pr_user) = &
1212                                    sums(k,pr_palm+1:pr_palm+max_pr_user) / &
1213                                    ngp_2dh_s_inner(k,sr)
1214          ENDDO
1215       ENDIF
1216
1217!
1218!--    Collect horizontal average in hom.
1219!--    Compute deduced averages (e.g. total heat flux)
1220       hom(:,1,3,sr)  = sums(:,3)      ! w
1221       hom(:,1,8,sr)  = sums(:,8)      ! e     profiles 5-7 are initial profiles
1222       hom(:,1,9,sr)  = sums(:,9)      ! km
1223       hom(:,1,10,sr) = sums(:,10)     ! kh
1224       hom(:,1,11,sr) = sums(:,11)     ! l
1225       hom(:,1,12,sr) = sums(:,12)     ! w"u"
1226       hom(:,1,13,sr) = sums(:,13)     ! w*u*
1227       hom(:,1,14,sr) = sums(:,14)     ! w"v"
1228       hom(:,1,15,sr) = sums(:,15)     ! w*v*
1229       hom(:,1,16,sr) = sums(:,16)     ! w"pt"
1230       hom(:,1,17,sr) = sums(:,17)     ! w*pt*
1231       hom(:,1,18,sr) = sums(:,16) + sums(:,17)    ! wpt
1232       hom(:,1,19,sr) = sums(:,12) + sums(:,13)    ! wu
1233       hom(:,1,20,sr) = sums(:,14) + sums(:,15)    ! wv
1234       hom(:,1,21,sr) = sums(:,21)     ! w*pt*BC
1235       hom(:,1,22,sr) = sums(:,16) + sums(:,21)    ! wptBC
1236                                       ! profile 24 is initial profile (sa)
1237                                       ! profiles 25-29 left empty for initial
1238                                       ! profiles
1239       hom(:,1,30,sr) = sums(:,30)     ! u*2
1240       hom(:,1,31,sr) = sums(:,31)     ! v*2
1241       hom(:,1,32,sr) = sums(:,32)     ! w*2
1242       hom(:,1,33,sr) = sums(:,33)     ! pt*2
1243       hom(:,1,34,sr) = sums(:,34)     ! e*
1244       hom(:,1,35,sr) = sums(:,35)     ! w*2pt*
1245       hom(:,1,36,sr) = sums(:,36)     ! w*pt*2
1246       hom(:,1,37,sr) = sums(:,37)     ! w*e*
1247       hom(:,1,38,sr) = sums(:,38)     ! w*3
1248       hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20_wp )**1.5_wp   ! Sw
1249       hom(:,1,40,sr) = sums(:,40)     ! p
1250       hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
1251       hom(:,1,46,sr) = sums(:,46)     ! w*vpt*       
1252       hom(:,1,47,sr) = sums(:,45) + sums(:,46)    ! wvpt
1253       hom(:,1,48,sr) = sums(:,48)     ! w"q" (w"qv")
1254       hom(:,1,49,sr) = sums(:,49)     ! w*q* (w*qv*)
1255       hom(:,1,50,sr) = sums(:,48) + sums(:,49)    ! wq (wqv)
1256       hom(:,1,51,sr) = sums(:,51)     ! w"qv"
1257       hom(:,1,52,sr) = sums(:,52)     ! w*qv*       
1258       hom(:,1,53,sr) = sums(:,52) + sums(:,51)    ! wq (wqv)
1259       hom(:,1,54,sr) = sums(:,54)     ! ql
1260       hom(:,1,55,sr) = sums(:,55)     ! w*u*u*/dz
1261       hom(:,1,56,sr) = sums(:,56)     ! w*p*/dz
1262       hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho )/dz
1263       hom(:,1,58,sr) = sums(:,58)     ! u"pt"
1264       hom(:,1,59,sr) = sums(:,59)     ! u*pt*
1265       hom(:,1,60,sr) = sums(:,58) + sums(:,59)    ! upt_t
1266       hom(:,1,61,sr) = sums(:,61)     ! v"pt"
1267       hom(:,1,62,sr) = sums(:,62)     ! v*pt*
1268       hom(:,1,63,sr) = sums(:,61) + sums(:,62)    ! vpt_t
1269       hom(:,1,64,sr) = sums(:,64)     ! rho
1270       hom(:,1,65,sr) = sums(:,65)     ! w"sa"
1271       hom(:,1,66,sr) = sums(:,66)     ! w*sa*
1272       hom(:,1,67,sr) = sums(:,65) + sums(:,66)    ! wsa
1273       hom(:,1,68,sr) = sums(:,68)     ! w*p*
1274       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho
1275       hom(:,1,70,sr) = sums(:,70)     ! q*2
1276       hom(:,1,71,sr) = sums(:,71)     ! prho
1277       hom(:,1,72,sr) = hyp * 1E-4_wp  ! hyp in dbar
1278       hom(:,1,73,sr) = sums(:,73)     ! nr
1279       hom(:,1,74,sr) = sums(:,74)     ! qr
1280       hom(:,1,75,sr) = sums(:,75)     ! qc
1281       hom(:,1,76,sr) = sums(:,76)     ! prr (precipitation rate)
1282                                       ! 77 is initial density profile
1283       hom(:,1,78,sr) = ug             ! ug
1284       hom(:,1,79,sr) = vg             ! vg
1285       hom(:,1,80,sr) = w_subs         ! w_subs
1286
1287       IF ( large_scale_forcing )  THEN
1288          hom(:,1,81,sr) = sums_ls_l(:,0)          ! td_lsa_lpt
1289          hom(:,1,82,sr) = sums_ls_l(:,1)          ! td_lsa_q
1290          IF ( use_subsidence_tendencies )  THEN
1291             hom(:,1,83,sr) = sums_ls_l(:,2)       ! td_sub_lpt
1292             hom(:,1,84,sr) = sums_ls_l(:,3)       ! td_sub_q
1293          ELSE
1294             hom(:,1,83,sr) = sums(:,83)           ! td_sub_lpt
1295             hom(:,1,84,sr) = sums(:,84)           ! td_sub_q
1296          ENDIF
1297          hom(:,1,85,sr) = sums(:,85)              ! td_nud_lpt
1298          hom(:,1,86,sr) = sums(:,86)              ! td_nud_q
1299          hom(:,1,87,sr) = sums(:,87)              ! td_nud_u
1300          hom(:,1,88,sr) = sums(:,88)              ! td_nud_v
1301       ENDIF
1302
1303       IF ( land_surface )  THEN
1304          hom(:,1,89,sr) = sums(:,89)              ! t_soil
1305                                                   ! 90 is initial t_soil profile
1306          hom(:,1,91,sr) = sums(:,91)              ! m_soil
1307                                                   ! 92 is initial m_soil profile
1308          hom(:,1,93,sr)  = sums(:,93)             ! ghf_eb
1309          hom(:,1,94,sr)  = sums(:,94)             ! shf_eb
1310          hom(:,1,95,sr)  = sums(:,95)             ! qsws_eb
1311          hom(:,1,96,sr)  = sums(:,96)             ! qsws_liq_eb
1312          hom(:,1,97,sr)  = sums(:,97)             ! qsws_soil_eb
1313          hom(:,1,98,sr)  = sums(:,98)             ! qsws_veg_eb
1314          hom(:,1,99,sr)  = sums(:,99)             ! r_a
1315          hom(:,1,100,sr) = sums(:,100)            ! r_s
1316
1317       ENDIF
1318
1319       IF ( radiation )  THEN
1320          hom(:,1,101 ,sr) = sums(:,101)            ! rad_net
1321          hom(:,1,102,sr)  = sums(:,102)            ! rad_sw_in
1322       ENDIF
1323
1324       hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
1325                                       ! upstream-parts u_x, u_y, u_z, v_x,
1326                                       ! v_y, usw. (in last but one profile)
1327       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
1328                                       ! u*, w'u', w'v', t* (in last profile)
1329
1330       IF ( max_pr_user > 0 )  THEN    ! user-defined profiles
1331          hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
1332                               sums(:,pr_palm+1:pr_palm+max_pr_user)
1333       ENDIF
1334
1335!
1336!--    Determine the boundary layer height using two different schemes.
1337!--    First scheme: Starting from the Earth's (Ocean's) surface, look for the
1338!--    first relative minimum (maximum) of the total heat flux.
1339!--    The corresponding height is assumed as the boundary layer height, if it
1340!--    is less than 1.5 times the height where the heat flux becomes negative
1341!--    (positive) for the first time.
1342       z_i(1) = 0.0_wp
1343       first = .TRUE.
1344
1345       IF ( ocean )  THEN
1346          DO  k = nzt, nzb+1, -1
1347             IF ( first .AND. hom(k,1,18,sr) < 0.0_wp                          &
1348                .AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp)  THEN
1349                first = .FALSE.
1350                height = zw(k)
1351             ENDIF
1352             IF ( hom(k,1,18,sr) < 0.0_wp  .AND.                               &
1353                  abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND.                        &
1354                  hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
1355                IF ( zw(k) < 1.5_wp * height )  THEN
1356                   z_i(1) = zw(k)
1357                ELSE
1358                   z_i(1) = height
1359                ENDIF
1360                EXIT
1361             ENDIF
1362          ENDDO
1363       ELSE
1364          DO  k = nzb, nzt-1
1365             IF ( first .AND. hom(k,1,18,sr) < 0.0_wp                          &
1366                .AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp )  THEN
1367                first = .FALSE.
1368                height = zw(k)
1369             ENDIF
1370             IF ( hom(k,1,18,sr) < 0.0_wp  .AND.                               & 
1371                  abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND.                        &
1372                  hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
1373                IF ( zw(k) < 1.5_wp * height )  THEN
1374                   z_i(1) = zw(k)
1375                ELSE
1376                   z_i(1) = height
1377                ENDIF
1378                EXIT
1379             ENDIF
1380          ENDDO
1381       ENDIF
1382
1383!
1384!--    Second scheme: Gradient scheme from Sullivan et al. (1998), modified
1385!--    by Uhlenbrock(2006). The boundary layer height is the height with the
1386!--    maximal local temperature gradient: starting from the second (the last
1387!--    but one) vertical gridpoint, the local gradient must be at least
1388!--    0.2K/100m and greater than the next four gradients.
1389!--    WARNING: The threshold value of 0.2K/100m must be adjusted for the
1390!--             ocean case!
1391       z_i(2) = 0.0_wp
1392       DO  k = nzb+1, nzt+1
1393          dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
1394       ENDDO
1395       dptdz_threshold = 0.2_wp / 100.0_wp
1396
1397       IF ( ocean )  THEN
1398          DO  k = nzt+1, nzb+5, -1
1399             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
1400                  dptdz(k) > dptdz(k-1)  .AND.  dptdz(k) > dptdz(k-2)  .AND.  &
1401                  dptdz(k) > dptdz(k-3)  .AND.  dptdz(k) > dptdz(k-4) )  THEN
1402                z_i(2) = zw(k-1)
1403                EXIT
1404             ENDIF
1405          ENDDO
1406       ELSE
1407          DO  k = nzb+1, nzt-3
1408             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
1409                  dptdz(k) > dptdz(k+1)  .AND.  dptdz(k) > dptdz(k+2)  .AND.  &
1410                  dptdz(k) > dptdz(k+3)  .AND.  dptdz(k) > dptdz(k+4) )  THEN
1411                z_i(2) = zw(k-1)
1412                EXIT
1413             ENDIF
1414          ENDDO
1415       ENDIF
1416
1417       hom(nzb+6,1,pr_palm,sr) = z_i(1)
1418       hom(nzb+7,1,pr_palm,sr) = z_i(2)
1419
1420!
1421!--    Computation of both the characteristic vertical velocity and
1422!--    the characteristic convective boundary layer temperature.
1423!--    The horizontal average at nzb+1 is input for the average temperature.
1424       IF ( hom(nzb,1,18,sr) > 0.0_wp .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8_wp  &
1425           .AND.  z_i(1) /= 0.0_wp )  THEN
1426          hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) *                 &
1427                                       hom(nzb,1,18,sr) *                      &
1428                                       ABS( z_i(1) ) )**0.333333333_wp
1429!--       so far this only works if Prandtl layer is used
1430          hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
1431       ELSE
1432          hom(nzb+8,1,pr_palm,sr)  = 0.0_wp
1433          hom(nzb+11,1,pr_palm,sr) = 0.0_wp
1434       ENDIF
1435
1436!
1437!--    Collect the time series quantities
1438       ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)     ! E
1439       ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)     ! E*
1440       ts_value(3,sr) = dt_3d
1441       ts_value(4,sr) = hom(nzb,1,pr_palm,sr)       ! u*
1442       ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr)     ! th*
1443       ts_value(6,sr) = u_max
1444       ts_value(7,sr) = v_max
1445       ts_value(8,sr) = w_max
1446       ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr)    ! new divergence
1447       ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr)    ! old Divergence
1448       ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr)    ! z_i(1)
1449       ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr)    ! z_i(2)
1450       ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr)    ! w*
1451       ts_value(14,sr) = hom(nzb,1,16,sr)           ! w'pt'   at k=0
1452       ts_value(15,sr) = hom(nzb+1,1,16,sr)         ! w'pt'   at k=1
1453       ts_value(16,sr) = hom(nzb+1,1,18,sr)         ! wpt     at k=1
1454       ts_value(17,sr) = hom(nzb,1,4,sr)            ! pt(0)
1455       ts_value(18,sr) = hom(nzb+1,1,4,sr)          ! pt(zp)
1456       ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr)    ! u'w'    at k=0
1457       ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)    ! v'w'    at k=0
1458       ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
1459
1460       IF ( ts_value(5,sr) /= 0.0_wp )  THEN
1461          ts_value(22,sr) = ts_value(4,sr)**2 / &
1462                            ( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
1463       ELSE
1464          ts_value(22,sr) = 10000.0_wp
1465       ENDIF
1466
1467       ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
1468
1469!
1470!--    Collect land surface model timeseries
1471       IF ( land_surface )  THEN
1472          ts_value(dots_soil  ,sr) = hom(nzb,1,93,sr)           ! ghf_eb
1473          ts_value(dots_soil+1,sr) = hom(nzb,1,94,sr)           ! shf_eb
1474          ts_value(dots_soil+2,sr) = hom(nzb,1,95,sr)           ! qsws_eb
1475          ts_value(dots_soil+3,sr) = hom(nzb,1,96,sr)           ! qsws_liq_eb
1476          ts_value(dots_soil+4,sr) = hom(nzb,1,97,sr)           ! qsws_soil_eb
1477          ts_value(dots_soil+5,sr) = hom(nzb,1,98,sr)           ! qsws_veg_eb
1478          ts_value(dots_soil+6,sr) = hom(nzb,1,99,sr)           ! r_a
1479          ts_value(dots_soil+7,sr) = hom(nzb,1,100,sr)          ! r_s
1480       ENDIF
1481!
1482!--    Collect radiation model timeseries
1483       IF ( radiation )  THEN
1484          ts_value(dots_rad,sr)   = hom(nzb,1,99,sr)           ! rad_net
1485          ts_value(dots_rad+1,sr) = hom(nzb,1,100,sr)          ! rad_sw_in
1486       ENDIF
1487
1488!
1489!--    Calculate additional statistics provided by the user interface
1490       CALL user_statistics( 'time_series', sr, 0 )
1491
1492    ENDDO    ! loop of the subregions
1493
1494!
1495!-- If required, sum up horizontal averages for subsequent time averaging
1496    IF ( do_sum )  THEN
1497       IF ( average_count_pr == 0 )  hom_sum = 0.0_wp
1498       hom_sum = hom_sum + hom(:,1,:,:)
1499       average_count_pr = average_count_pr + 1
1500       do_sum = .FALSE.
1501    ENDIF
1502
1503!
1504!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
1505!-- This flag is reset after each time step in time_integration.
1506    flow_statistics_called = .TRUE.
1507
1508    CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
1509
1510
1511 END SUBROUTINE flow_statistics
1512
1513
1514#else
1515
1516
1517!------------------------------------------------------------------------------!
1518! flow statistics - accelerator version
1519!------------------------------------------------------------------------------!
1520 SUBROUTINE flow_statistics
1521
1522    USE arrays_3d,                                                             &
1523        ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, p, prho, pt, q, qc, ql, qr, qs, &
1524               qsws, qswst, rho, sa, saswsb, saswst, shf, td_lsa_lpt, td_lsa_q,&
1525               td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us, usws,    &
1526               uswst, vsws, v, vg, vpt, vswst, w, w_subs, zw
1527                 
1528       
1529    USE cloud_parameters,                                                      &
1530        ONLY:  l_d_cp, prr, pt_d_t
1531       
1532    USE control_parameters,                                                    &
1533        ONLY :  average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
1534                dt_3d, g, humidity, icloud_scheme, kappa, large_scale_forcing, &
1535                large_scale_subsidence, max_pr_user, message_string,           &
1536                monotonic_adjustment, 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. .NOT. monotonic_adjustment .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. monotonic_adjustment .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. monotonic_adjustment .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. monotonic_adjustment &
2672                                     .OR.  sr /= 0 ) )  THEN
2673
2674          !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l, w ) create( s1 )
2675          DO  k = nzb, nzt_diff
2676             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2677             DO  i = nxl, nxr
2678                DO  j = nys, nyn
2679                   s1 = s1 + 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +          &
2680                                        q(k+1,j,i) - hom(k+1,1,41,sr) )        &
2681                                    * w(k,j,i) * rmask(j,i,sr)                 &
2682                                    * rflags_invers(j,i,k+1)
2683                ENDDO
2684             ENDDO
2685             sums_l(k,49,tn) = s1
2686          ENDDO
2687          !$acc end parallel loop
2688
2689       ENDIF
2690
2691!
2692!--    For speed optimization fluxes which have been computed in part directly
2693!--    inside the WS advection routines are treated seperatly
2694!--    Momentum fluxes first:
2695       IF ( .NOT. ws_scheme_mom  .OR.  sr /= 0  )  THEN
2696
2697          !$OMP DO
2698          !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, u, v, w ) create( s1, s2 )
2699          DO  k = nzb, nzt_diff
2700             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2701             DO  i = nxl, nxr
2702                DO  j = nys, nyn
2703                   ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +               &
2704                                    u(k+1,j,i) - hom(k+1,1,1,sr) )
2705                   vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +               &
2706                                    v(k+1,j,i) - hom(k+1,1,2,sr) )
2707!
2708!--                Momentum flux w*u*
2709                   s1 = s1 + 0.5_wp * ( w(k,j,i-1) + w(k,j,i) )                &
2710                                    * ust * rmask(j,i,sr)                      &
2711                                    * rflags_invers(j,i,k+1)
2712!
2713!--                Momentum flux w*v*
2714                   s2 = s2 + 0.5_wp * ( w(k,j-1,i) + w(k,j,i) )                &
2715                                    * vst * rmask(j,i,sr)                      & 
2716                                    * rflags_invers(j,i,k+1)
2717                ENDDO
2718             ENDDO
2719             sums_l(k,13,tn) = s1
2720             sums_l(k,15,tn) = s1
2721          ENDDO
2722          !$acc end parallel loop
2723
2724       ENDIF
2725
2726       IF ( .NOT. ws_scheme_sca .OR. monotonic_adjustment .OR.  sr /= 0 )  THEN
2727
2728          !$OMP DO
2729          !$acc parallel loop gang present( hom, pt, rflags_invers, rmask, sums_l, w ) create( s1 )
2730          DO  k = nzb, nzt_diff
2731             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2732             DO  i = nxl, nxr
2733                DO  j = nys, nyn
2734!
2735!--                Vertical heat flux
2736                   s1 = s1 + 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) +          &
2737                                        pt(k+1,j,i) - hom(k+1,1,4,sr) )        &
2738                                    * w(k,j,i) * rmask(j,i,sr)                 & 
2739                                    * rflags_invers(j,i,k+1)
2740                ENDDO
2741             ENDDO
2742             sums_l(k,17,tn) = s1
2743          ENDDO
2744          !$acc end parallel loop
2745
2746          IF ( humidity )  THEN
2747
2748             !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l, w ) create( s1 )
2749             DO  k = nzb, nzt_diff
2750                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2751                DO  i = nxl, nxr
2752                   DO  j = nys, nyn
2753                      s1 = s1 + 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +       &
2754                                           q(k+1,j,i) - hom(k+1,1,41,sr) )     &
2755                                       * w(k,j,i) * rmask(j,i,sr)              &
2756                                       * rflags_invers(j,i,k+1)
2757                   ENDDO
2758                ENDDO
2759                sums_l(k,49,tn) = s1
2760             ENDDO
2761             !$acc end parallel loop
2762
2763          ENDIF
2764
2765       ENDIF
2766
2767
2768!
2769!--    Density at top follows Neumann condition
2770       IF ( ocean )  THEN
2771          !$acc parallel present( sums_l )
2772          sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
2773          sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
2774          !$acc end parallel
2775       ENDIF
2776
2777!
2778!--    Divergence of vertical flux of resolved scale energy and pressure
2779!--    fluctuations as well as flux of pressure fluctuation itself (68).
2780!--    First calculate the products, then the divergence.
2781!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
2782       IF ( hom(nzb+1,2,55,0) /= 0.0_wp  .OR.  hom(nzb+1,2,68,0) /= 0.0_wp )  THEN
2783
2784          STOP '+++ openACC porting for vertical flux div of resolved scale TKE in flow_statistics is still missing'
2785          sums_ll = 0.0_wp  ! local array
2786
2787          !$OMP DO
2788          DO  i = nxl, nxr
2789             DO  j = nys, nyn
2790                DO  k = nzb_s_inner(j,i)+1, nzt
2791
2792                   sums_ll(k,1) = sums_ll(k,1) + 0.5_wp * w(k,j,i) * (         &
2793                  ( 0.25_wp * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1)    &
2794                              - 0.5_wp * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) )   &
2795                              ) )**2                                           &
2796                + ( 0.25_wp * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i)    &
2797                              - 0.5_wp * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) )   &
2798                              ) )**2                                           &
2799                + w(k,j,i)**2                                        )
2800
2801                   sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i)             &
2802                                               * ( p(k,j,i) + p(k+1,j,i) )
2803
2804                ENDDO
2805             ENDDO
2806          ENDDO
2807          sums_ll(0,1)     = 0.0_wp    ! because w is zero at the bottom
2808          sums_ll(nzt+1,1) = 0.0_wp
2809          sums_ll(0,2)     = 0.0_wp
2810          sums_ll(nzt+1,2) = 0.0_wp
2811
2812          DO  k = nzb+1, nzt
2813             sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
2814             sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
2815             sums_l(k,68,tn) = sums_ll(k,2)
2816          ENDDO
2817          sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
2818          sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
2819          sums_l(nzb,68,tn) = 0.0_wp    ! because w* = 0 at nzb
2820
2821       ENDIF
2822
2823!
2824!--    Divergence of vertical flux of SGS TKE and the flux itself (69)
2825       IF ( hom(nzb+1,2,57,0) /= 0.0_wp  .OR.  hom(nzb+1,2,69,0) /= 0.0_wp )  THEN
2826
2827          STOP '+++ openACC porting for vertical flux div of SGS TKE in flow_statistics is still missing'
2828          !$OMP DO
2829          DO  i = nxl, nxr
2830             DO  j = nys, nyn
2831                DO  k = nzb_s_inner(j,i)+1, nzt
2832
2833                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5_wp * (              &
2834                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
2835                 - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
2836                                                                ) * ddzw(k)
2837
2838                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5_wp * (              &
2839                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
2840                                                                )
2841
2842                ENDDO
2843             ENDDO
2844          ENDDO
2845          sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
2846          sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
2847
2848       ENDIF
2849
2850!
2851!--    Horizontal heat fluxes (subgrid, resolved, total).
2852!--    Do it only, if profiles shall be plotted.
2853       IF ( hom(nzb+1,2,58,0) /= 0.0_wp ) THEN
2854
2855          STOP '+++ openACC porting for horizontal flux calculation in flow_statistics is still missing'
2856          !$OMP DO
2857          DO  i = nxl, nxr
2858             DO  j = nys, nyn
2859                DO  k = nzb_s_inner(j,i)+1, nzt
2860!
2861!--                Subgrid horizontal heat fluxes u"pt", v"pt"
2862                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5_wp *                &
2863                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
2864                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
2865                                                 * ddx * rmask(j,i,sr)
2866                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
2867                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
2868                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
2869                                                 * ddy * rmask(j,i,sr)
2870!
2871!--                Resolved horizontal heat fluxes u*pt*, v*pt*
2872                   sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
2873                                     ( u(k,j,i) - hom(k,1,1,sr) ) * 0.5_wp *   &
2874                                     ( pt(k,j,i-1) - hom(k,1,4,sr) +           &
2875                                       pt(k,j,i)   - hom(k,1,4,sr) )
2876                   pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
2877                                    pt(k,j,i)   - hom(k,1,4,sr) )
2878                   sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
2879                                     ( v(k,j,i) - hom(k,1,2,sr) ) * 0.5_wp *   &
2880                                     ( pt(k,j-1,i) - hom(k,1,4,sr) +           &
2881                                       pt(k,j,i)   - hom(k,1,4,sr) )
2882                ENDDO
2883             ENDDO
2884          ENDDO
2885!
2886!--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
2887          sums_l(nzb,58,tn) = 0.0_wp
2888          sums_l(nzb,59,tn) = 0.0_wp
2889          sums_l(nzb,60,tn) = 0.0_wp
2890          sums_l(nzb,61,tn) = 0.0_wp
2891          sums_l(nzb,62,tn) = 0.0_wp
2892          sums_l(nzb,63,tn) = 0.0_wp
2893
2894       ENDIF
2895
2896!
2897!--    Collect current large scale advection and subsidence tendencies for
2898!--    data output
2899       IF ( large_scale_forcing  .AND.  ( simulated_time .GT. 0.0_wp ) )  THEN
2900!
2901!--       Interpolation in time of LSF_DATA
2902          nt = 1
2903          DO WHILE ( simulated_time - dt_3d > time_vert(nt) )
2904             nt = nt + 1
2905          ENDDO
2906          IF ( simulated_time - dt_3d /= time_vert(nt) )  THEN
2907            nt = nt - 1
2908          ENDIF
2909
2910          fac = ( simulated_time - dt_3d - time_vert(nt) )                     &
2911                / ( time_vert(nt+1)-time_vert(nt) )
2912
2913
2914          DO  k = nzb, nzt
2915             sums_ls_l(k,0) = td_lsa_lpt(k,nt)                                 &
2916                              + fac * ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )
2917             sums_ls_l(k,1) = td_lsa_q(k,nt)                                   &
2918                              + fac * ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) )
2919          ENDDO
2920
2921          sums_ls_l(nzt+1,0) = sums_ls_l(nzt,0)
2922          sums_ls_l(nzt+1,1) = sums_ls_l(nzt,1)
2923
2924          IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
2925
2926             DO  k = nzb, nzt
2927                sums_ls_l(k,2) = td_sub_lpt(k,nt) + fac *                      &
2928                                 ( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )
2929                sums_ls_l(k,3) = td_sub_q(k,nt) + fac *                        &
2930                                 ( td_sub_q(k,nt+1) - td_sub_q(k,nt) )
2931             ENDDO
2932
2933             sums_ls_l(nzt+1,2) = sums_ls_l(nzt,2)
2934             sums_ls_l(nzt+1,3) = sums_ls_l(nzt,3)
2935
2936          ENDIF
2937
2938       ENDIF
2939
2940
2941       IF ( land_surface )  THEN
2942          !$OMP DO
2943          DO  i = nxl, nxr
2944             DO  j =  nys, nyn
2945                DO  k = nzb_soil, nzt_soil
2946                   sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil(k,j,i) * rmask(j,i,sr)
2947                   sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil(k,j,i) * rmask(j,i,sr)
2948                ENDDO
2949             ENDDO
2950          ENDDO
2951       ENDIF
2952
2953
2954!
2955!--    Calculate the user-defined profiles
2956       CALL user_statistics( 'profiles', sr, tn )
2957       !$OMP END PARALLEL
2958
2959!
2960!--    Summation of thread sums
2961       IF ( threads_per_task > 1 )  THEN
2962          STOP '+++ openACC porting for threads_per_task > 1 in flow_statistics is still missing'
2963          DO  i = 1, threads_per_task-1
2964             sums_l(:,3,0)          = sums_l(:,3,0) + sums_l(:,3,i)
2965             sums_l(:,4:40,0)       = sums_l(:,4:40,0) + sums_l(:,4:40,i)
2966             sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
2967                                      sums_l(:,45:pr_palm,i)
2968             IF ( max_pr_user > 0 )  THEN
2969                sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
2970                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
2971                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
2972             ENDIF
2973          ENDDO
2974       ENDIF
2975
2976       !$acc update host( hom, sums, sums_l )
2977
2978#if defined( __parallel )
2979
2980!
2981!--    Compute total sum from local sums
2982       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2983       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, &
2984                           MPI_SUM, comm2d, ierr )
2985       IF ( large_scale_forcing )  THEN
2986          CALL MPI_ALLREDUCE( sums_ls_l(nzb,2), sums(nzb,83), ngp_sums_ls,     &
2987                              MPI_REAL, MPI_SUM, comm2d, ierr )
2988       ENDIF
2989#else
2990       sums = sums_l(:,:,0)
2991       IF ( large_scale_forcing )  THEN
2992          sums(:,81:88) = sums_ls_l
2993       ENDIF
2994#endif
2995
2996!
2997!--    Final values are obtained by division by the total number of grid points
2998!--    used for summation. After that store profiles.
2999!--    Profiles:
3000       DO  k = nzb, nzt+1
3001          sums(k,3)               = sums(k,3)           / ngp_2dh(sr)
3002          sums(k,8:11)            = sums(k,8:11)        / ngp_2dh_s_inner(k,sr)
3003          sums(k,12:22)           = sums(k,12:22)       / ngp_2dh(sr)
3004          sums(k,23:29)           = sums(k,23:29)       / ngp_2dh_s_inner(k,sr)
3005          sums(k,30:32)           = sums(k,30:32)       / ngp_2dh(sr)
3006          sums(k,33:34)           = sums(k,33:34)       / ngp_2dh_s_inner(k,sr)
3007          sums(k,35:39)           = sums(k,35:39)       / ngp_2dh(sr)
3008          sums(k,40)              = sums(k,40)          / ngp_2dh_s_inner(k,sr)
3009          sums(k,45:53)           = sums(k,45:53)       / ngp_2dh(sr)
3010          sums(k,54)              = sums(k,54)          / ngp_2dh_s_inner(k,sr)
3011          sums(k,55:63)           = sums(k,55:63)       / ngp_2dh(sr)
3012          sums(k,64)              = sums(k,64)          / ngp_2dh_s_inner(k,sr)
3013          sums(k,70:80)           = sums(k,70:80)       / ngp_2dh_s_inner(k,sr)
3014          sums(k,81:88)           = sums(k,81:88)       / ngp_2dh(sr)
3015          sums(k,89:100)           = sums(k,89:100)     / ngp_2dh(sr)
3016          sums(k,101:pr_palm-2)    = sums(k,101:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
3017       ENDDO
3018
3019!--    Upstream-parts
3020       sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
3021!--    u* and so on
3022!--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
3023!--    size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
3024!--    above the topography, they are being divided by ngp_2dh(sr)
3025       sums(nzb:nzb+3,pr_palm)    = sums(nzb:nzb+3,pr_palm)    / &
3026                                    ngp_2dh(sr)
3027       sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)       / &    ! qs
3028                                    ngp_2dh(sr)
3029!--    eges, e*
3030       sums(nzb+4:nzb+5,pr_palm)  = sums(nzb+4:nzb+5,pr_palm)  / &
3031                                    ngp_3d(sr)
3032!--    Old and new divergence
3033       sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
3034                                    ngp_3d_inner(sr)
3035
3036!--    User-defined profiles
3037       IF ( max_pr_user > 0 )  THEN
3038          DO  k = nzb, nzt+1
3039             sums(k,pr_palm+1:pr_palm+max_pr_user) = &
3040                                    sums(k,pr_palm+1:pr_palm+max_pr_user) / &
3041                                    ngp_2dh_s_inner(k,sr)
3042          ENDDO
3043       ENDIF
3044
3045!
3046!--    Collect horizontal average in hom.
3047!--    Compute deduced averages (e.g. total heat flux)
3048       hom(:,1,3,sr)  = sums(:,3)      ! w
3049       hom(:,1,8,sr)  = sums(:,8)      ! e     profiles 5-7 are initial profiles
3050       hom(:,1,9,sr)  = sums(:,9)      ! km
3051       hom(:,1,10,sr) = sums(:,10)     ! kh
3052       hom(:,1,11,sr) = sums(:,11)     ! l
3053       hom(:,1,12,sr) = sums(:,12)     ! w"u"
3054       hom(:,1,13,sr) = sums(:,13)     ! w*u*
3055       hom(:,1,14,sr) = sums(:,14)     ! w"v"
3056       hom(:,1,15,sr) = sums(:,15)     ! w*v*
3057       hom(:,1,16,sr) = sums(:,16)     ! w"pt"
3058       hom(:,1,17,sr) = sums(:,17)     ! w*pt*
3059       hom(:,1,18,sr) = sums(:,16) + sums(:,17)    ! wpt
3060       hom(:,1,19,sr) = sums(:,12) + sums(:,13)    ! wu
3061       hom(:,1,20,sr) = sums(:,14) + sums(:,15)    ! wv
3062       hom(:,1,21,sr) = sums(:,21)     ! w*pt*BC
3063       hom(:,1,22,sr) = sums(:,16) + sums(:,21)    ! wptBC
3064                                       ! profile 24 is initial profile (sa)
3065                                       ! profiles 25-29 left empty for initial
3066                                       ! profiles
3067       hom(:,1,30,sr) = sums(:,30)     ! u*2
3068       hom(:,1,31,sr) = sums(:,31)     ! v*2
3069       hom(:,1,32,sr) = sums(:,32)     ! w*2
3070       hom(:,1,33,sr) = sums(:,33)     ! pt*2
3071       hom(:,1,34,sr) = sums(:,34)     ! e*
3072       hom(:,1,35,sr) = sums(:,35)     ! w*2pt*
3073       hom(:,1,36,sr) = sums(:,36)     ! w*pt*2
3074       hom(:,1,37,sr) = sums(:,37)     ! w*e*
3075       hom(:,1,38,sr) = sums(:,38)     ! w*3
3076       hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20_wp )**1.5_wp   ! Sw
3077       hom(:,1,40,sr) = sums(:,40)     ! p
3078       hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
3079       hom(:,1,46,sr) = sums(:,46)     ! w*vpt*
3080       hom(:,1,47,sr) = sums(:,45) + sums(:,46)    ! wvpt
3081       hom(:,1,48,sr) = sums(:,48)     ! w"q" (w"qv")
3082       hom(:,1,49,sr) = sums(:,49)     ! w*q* (w*qv*)
3083       hom(:,1,50,sr) = sums(:,48) + sums(:,49)    ! wq (wqv)
3084       hom(:,1,51,sr) = sums(:,51)     ! w"qv"
3085       hom(:,1,52,sr) = sums(:,52)     ! w*qv*
3086       hom(:,1,53,sr) = sums(:,52) + sums(:,51)    ! wq (wqv)
3087       hom(:,1,54,sr) = sums(:,54)     ! ql
3088       hom(:,1,55,sr) = sums(:,55)     ! w*u*u*/dz
3089       hom(:,1,56,sr) = sums(:,56)     ! w*p*/dz
3090       hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho )/dz
3091       hom(:,1,58,sr) = sums(:,58)     ! u"pt"
3092       hom(:,1,59,sr) = sums(:,59)     ! u*pt*
3093       hom(:,1,60,sr) = sums(:,58) + sums(:,59)    ! upt_t
3094       hom(:,1,61,sr) = sums(:,61)     ! v"pt"
3095       hom(:,1,62,sr) = sums(:,62)     ! v*pt*
3096       hom(:,1,63,sr) = sums(:,61) + sums(:,62)    ! vpt_t
3097       hom(:,1,64,sr) = sums(:,64)     ! rho
3098       hom(:,1,65,sr) = sums(:,65)     ! w"sa"
3099       hom(:,1,66,sr) = sums(:,66)     ! w*sa*
3100       hom(:,1,67,sr) = sums(:,65) + sums(:,66)    ! wsa
3101       hom(:,1,68,sr) = sums(:,68)     ! w*p*
3102       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho
3103       hom(:,1,70,sr) = sums(:,70)     ! q*2
3104       hom(:,1,71,sr) = sums(:,71)     ! prho
3105       hom(:,1,72,sr) = hyp * 1E-4_wp     ! hyp in dbar
3106       hom(:,1,73,sr) = sums(:,73)     ! nr
3107       hom(:,1,74,sr) = sums(:,74)     ! qr
3108       hom(:,1,75,sr) = sums(:,75)     ! qc
3109       hom(:,1,76,sr) = sums(:,76)     ! prr (precipitation rate)
3110                                       ! 77 is initial density profile
3111       hom(:,1,78,sr) = ug             ! ug
3112       hom(:,1,79,sr) = vg             ! vg
3113       hom(:,1,80,sr) = w_subs         ! w_subs
3114
3115       IF ( large_scale_forcing )  THEN
3116          hom(:,1,81,sr) = sums_ls_l(:,0)          ! td_lsa_lpt
3117          hom(:,1,82,sr) = sums_ls_l(:,1)          ! td_lsa_q
3118          IF ( use_subsidence_tendencies )  THEN
3119             hom(:,1,83,sr) = sums_ls_l(:,2)       ! td_sub_lpt
3120             hom(:,1,84,sr) = sums_ls_l(:,3)       ! td_sub_q
3121          ELSE
3122             hom(:,1,83,sr) = sums(:,83)           ! td_sub_lpt
3123             hom(:,1,84,sr) = sums(:,84)           ! td_sub_q
3124          ENDIF
3125          hom(:,1,85,sr) = sums(:,85)              ! td_nud_lpt
3126          hom(:,1,86,sr) = sums(:,86)              ! td_nud_q
3127          hom(:,1,87,sr) = sums(:,87)              ! td_nud_u
3128          hom(:,1,88,sr) = sums(:,88)              ! td_nud_v
3129       END IF
3130
3131       hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
3132                                       ! upstream-parts u_x, u_y, u_z, v_x,
3133                                       ! v_y, usw. (in last but one profile)
3134       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
3135                                       ! u*, w'u', w'v', t* (in last profile)
3136
3137       IF ( max_pr_user > 0 )  THEN    ! user-defined profiles
3138          hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
3139                               sums(:,pr_palm+1:pr_palm+max_pr_user)
3140       ENDIF
3141
3142!
3143!--    Determine the boundary layer height using two different schemes.
3144!--    First scheme: Starting from the Earth's (Ocean's) surface, look for the
3145!--    first relative minimum (maximum) of the total heat flux.
3146!--    The corresponding height is assumed as the boundary layer height, if it
3147!--    is less than 1.5 times the height where the heat flux becomes negative
3148!--    (positive) for the first time.
3149       z_i(1) = 0.0_wp
3150       first = .TRUE.
3151
3152       IF ( ocean )  THEN
3153          DO  k = nzt, nzb+1, -1
3154             IF ( first .AND. hom(k,1,18,sr) < 0.0_wp                          &
3155                .AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp )  THEN
3156                first = .FALSE.
3157                height = zw(k)
3158             ENDIF
3159             IF ( hom(k,1,18,sr) < 0.0_wp  .AND.                               &
3160                  abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND.                        &
3161                  hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
3162                IF ( zw(k) < 1.5_wp * height )  THEN
3163                   z_i(1) = zw(k)
3164                ELSE
3165                   z_i(1) = height
3166                ENDIF
3167                EXIT
3168             ENDIF
3169          ENDDO
3170       ELSE
3171          DO  k = nzb, nzt-1
3172             IF ( first .AND. hom(k,1,18,sr) < 0.0_wp                          &
3173                .AND. abs(hom(k,1,18,sr)) > 1.0E-8_wp )  THEN
3174                first = .FALSE.
3175                height = zw(k)
3176             ENDIF
3177             IF ( hom(k,1,18,sr) < 0.0  .AND. &
3178                  abs(hom(k,1,18,sr)) > 1.0E-8_wp .AND. &
3179                  hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
3180                IF ( zw(k) < 1.5_wp * height )  THEN
3181                   z_i(1) = zw(k)
3182                ELSE
3183                   z_i(1) = height
3184                ENDIF
3185                EXIT
3186             ENDIF
3187          ENDDO
3188       ENDIF
3189
3190!
3191!--    Second scheme: Gradient scheme from Sullivan et al. (1998), modified
3192!--    by Uhlenbrock(2006). The boundary layer height is the height with the
3193!--    maximal local temperature gradient: starting from the second (the last
3194!--    but one) vertical gridpoint, the local gradient must be at least
3195!--    0.2K/100m and greater than the next four gradients.
3196!--    WARNING: The threshold value of 0.2K/100m must be adjusted for the
3197!--             ocean case!
3198       z_i(2) = 0.0_wp
3199       DO  k = nzb+1, nzt+1
3200          dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
3201       ENDDO
3202       dptdz_threshold = 0.2_wp / 100.0_wp
3203
3204       IF ( ocean )  THEN
3205          DO  k = nzt+1, nzb+5, -1
3206             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
3207                  dptdz(k) > dptdz(k-1)  .AND.  dptdz(k) > dptdz(k-2)  .AND.  &
3208                  dptdz(k) > dptdz(k-3)  .AND.  dptdz(k) > dptdz(k-4) )  THEN
3209                z_i(2) = zw(k-1)
3210                EXIT
3211             ENDIF
3212          ENDDO
3213       ELSE
3214          DO  k = nzb+1, nzt-3
3215             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
3216                  dptdz(k) > dptdz(k+1)  .AND.  dptdz(k) > dptdz(k+2)  .AND.  &
3217                  dptdz(k) > dptdz(k+3)  .AND.  dptdz(k) > dptdz(k+4) )  THEN
3218                z_i(2) = zw(k-1)
3219                EXIT
3220             ENDIF
3221          ENDDO
3222       ENDIF
3223
3224       hom(nzb+6,1,pr_palm,sr) = z_i(1)
3225       hom(nzb+7,1,pr_palm,sr) = z_i(2)
3226
3227!
3228!--    Computation of both the characteristic vertical velocity and
3229!--    the characteristic convective boundary layer temperature.
3230!--    The horizontal average at nzb+1 is input for the average temperature.
3231       IF ( hom(nzb,1,18,sr) > 0.0_wp .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8_wp  &
3232           .AND.  z_i(1) /= 0.0_wp )  THEN
3233          hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) *                 &
3234                                       hom(nzb,1,18,sr) *                      &
3235                                       ABS( z_i(1) ) )**0.333333333_wp
3236!--       so far this only works if Prandtl layer is used
3237          hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
3238       ELSE
3239          hom(nzb+8,1,pr_palm,sr)  = 0.0_wp
3240          hom(nzb+11,1,pr_palm,sr) = 0.0_wp
3241       ENDIF
3242
3243!
3244!--    Collect the time series quantities
3245       ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)     ! E
3246       ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)     ! E*
3247       ts_value(3,sr) = dt_3d
3248       ts_value(4,sr) = hom(nzb,1,pr_palm,sr)       ! u*
3249       ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr)     ! th*
3250       ts_value(6,sr) = u_max
3251       ts_value(7,sr) = v_max
3252       ts_value(8,sr) = w_max
3253       ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr)    ! new divergence
3254       ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr)    ! old Divergence
3255       ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr)    ! z_i(1)
3256       ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr)    ! z_i(2)
3257       ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr)    ! w*
3258       ts_value(14,sr) = hom(nzb,1,16,sr)           ! w'pt'   at k=0
3259       ts_value(15,sr) = hom(nzb+1,1,16,sr)         ! w'pt'   at k=1
3260       ts_value(16,sr) = hom(nzb+1,1,18,sr)         ! wpt     at k=1
3261       ts_value(17,sr) = hom(nzb,1,4,sr)            ! pt(0)
3262       ts_value(18,sr) = hom(nzb+1,1,4,sr)          ! pt(zp)
3263       ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr)    ! u'w'    at k=0
3264       ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)    ! v'w'    at k=0
3265       ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
3266
3267       IF ( ts_value(5,sr) /= 0.0_wp )  THEN
3268          ts_value(22,sr) = ts_value(4,sr)**2_wp / &
3269                            ( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
3270       ELSE
3271          ts_value(22,sr) = 10000.0_wp
3272       ENDIF
3273
3274       ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
3275
3276!
3277!--    Collect land surface model timeseries
3278       IF ( land_surface )  THEN
3279          ts_value(dots_soil  ,sr) = hom(nzb,1,93,sr)           ! ghf_eb
3280          ts_value(dots_soil+1,sr) = hom(nzb,1,94,sr)           ! shf_eb
3281          ts_value(dots_soil+2,sr) = hom(nzb,1,95,sr)           ! qsws_eb
3282          ts_value(dots_soil+3,sr) = hom(nzb,1,96,sr)           ! qsws_liq_eb
3283          ts_value(dots_soil+4,sr) = hom(nzb,1,97,sr)           ! qsws_soil_eb
3284          ts_value(dots_soil+5,sr) = hom(nzb,1,98,sr)           ! qsws_veg_eb
3285          ts_value(dots_soil+6,sr) = hom(nzb,1,99,sr)           ! r_a
3286          ts_value(dots_soil+7,sr) = hom(nzb,1,100,sr)          ! r_s
3287       ENDIF
3288!
3289!--    Collect radiation model timeseries
3290       IF ( radiation )  THEN
3291          ts_value(dots_rad,sr)   = hom(nzb,1,99,sr)           ! rad_net
3292          ts_value(dots_rad+1,sr) = hom(nzb,1,100,sr)          ! rad_sw_in
3293       ENDIF
3294
3295!
3296!--    Calculate additional statistics provided by the user interface
3297       CALL user_statistics( 'time_series', sr, 0 )
3298
3299    ENDDO    ! loop of the subregions
3300
3301    !$acc end data
3302
3303!
3304!-- If required, sum up horizontal averages for subsequent time averaging
3305    IF ( do_sum )  THEN
3306       IF ( average_count_pr == 0 )  hom_sum = 0.0_wp
3307       hom_sum = hom_sum + hom(:,1,:,:)
3308       average_count_pr = average_count_pr + 1
3309       do_sum = .FALSE.
3310    ENDIF
3311
3312!
3313!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
3314!-- This flag is reset after each time step in time_integration.
3315    flow_statistics_called = .TRUE.
3316
3317    CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
3318
3319
3320 END SUBROUTINE flow_statistics
3321#endif
Note: See TracBrowser for help on using the repository browser.