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

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

bugfix: "copyin" replaced by "update device" in openacc-branch

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