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

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

minor changes in profile data output of lsf tendencies, variables renamed

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