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

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

further modularization of land surface model (2D/3D output and restart data). Bugfix for restart runs without land surface model

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