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

Last change on this file since 1322 was 1322, checked in by raasch, 7 years ago

REAL functions and a lot of REAL constants provided with KIND-attribute,
some small bugfixes

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