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

Last change on this file since 1652 was 1652, checked in by raasch, 6 years ago

bugfix in calculation of vertical profile of energy production by turbulent transport of TKE
bugfix for output of configuration settings (mbuild)
output of login init command (mbuild)

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