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

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

openacc adjustments and bugfixes

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