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

Last change on this file since 1450 was 1450, checked in by heinze, 7 years ago

bugfixes for large_scale forcing and nudging

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