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

Last change on this file since 2032 was 2032, checked in by knoop, 5 years ago

last commit documented

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