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

Last change on this file since 1585 was 1585, checked in by maronga, 6 years ago

Added support for RRTMG radiation code

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