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

Last change on this file since 1387 was 1387, checked in by boeske, 7 years ago

last commit documented

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