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

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

last commit documented

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