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

Last change on this file since 1323 was 1323, checked in by raasch, 10 years ago

last commit documented

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