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

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

several bugfixes related to the new surface layer routine and land-surface-radiation interaction

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