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

Last change on this file since 1654 was 1654, checked in by raasch, 6 years ago

FORTRAN bugfix of r1652

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