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

Last change on this file since 1960 was 1960, checked in by suehring, 5 years ago

Separate balance equations for humidity and passive_scalar

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