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

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

last commit documented / copyright update

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