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

Last change on this file since 1366 was 1366, checked in by boeske, 10 years ago

last commit documented

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