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

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

last commit documented

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