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

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

Forced header and separation lines into 80 columns

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