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

Last change on this file since 1822 was 1822, checked in by hoffmann, 5 years ago

changes in LPM and bulk cloud microphysics

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