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

Last change on this file since 1498 was 1498, checked in by suehring, 7 years ago

Comments added

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