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

Last change on this file since 1314 was 1310, checked in by raasch, 11 years ago

update of GPL copyright

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