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

Last change on this file since 1851 was 1851, checked in by maronga, 5 years ago

last commit documented

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