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

Last change on this file since 1683 was 1683, checked in by knoop, 6 years ago

last commit documented

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