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

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

Anelastic approximation implemented

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