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

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

Last commit documented

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