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

Last change on this file since 1483 was 1483, checked in by raasch, 7 years ago

last commit documented

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