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

Last change on this file since 1748 was 1748, checked in by raasch, 5 years ago

last commit documented

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