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

Last change on this file since 2001 was 2001, checked in by knoop, 5 years ago

last commit documented

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