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

Last change on this file since 1738 was 1738, checked in by raasch, 5 years ago

bugfixes for calculations in statistical regions which do not contain grid points in the lowest vertical levels, mean surface level height considered in the calculation of the characteristic vertical velocity

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