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

Last change on this file since 1669 was 1659, checked in by raasch, 9 years ago

last commit documented

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