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

Last change on this file since 1556 was 1556, checked in by maronga, 6 years ago

last commit documented

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