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

Last change on this file since 4338 was 4329, checked in by motisi, 5 years ago

Renamed wall_flags_0 to wall_flags_static_0

  • Property svn:keywords set to Id
File size: 106.1 KB
Line 
1!> @file flow_statistics.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2019 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: flow_statistics.f90 4329 2019-12-10 15:46:36Z suehring $
27! Renamed wall_flags_0 to wall_flags_static_0
28!
29! 4182 2019-08-22 15:20:23Z scharf
30! Corrected "Former revisions" section
31!
32! 4131 2019-08-02 11:06:18Z monakurppa
33! Allow profile output for salsa variables.
34!
35! 4039 2019-06-18 10:32:41Z suehring
36! Correct conversion to kinematic scalar fluxes in case of pw-scheme and
37! statistic regions
38!
39! 3828 2019-03-27 19:36:23Z raasch
40! unused variables removed
41!
42! 3676 2019-01-16 15:07:05Z knoop
43! Bugfix, terminate OMP Parallel block
44!
45! Revision 1.1  1997/08/11 06:15:17  raasch
46! Initial revision
47!
48!
49! Description:
50! ------------
51!> Compute average profiles and further average flow quantities for the different
52!> user-defined (sub-)regions. The region indexed 0 is the total model domain.
53!>
54!> @note For simplicity, nzb_s_inner and nzb_diff_s_inner are being used as a
55!>       lower vertical index for k-loops for all variables, although strictly
56!>       speaking the k-loops would have to be split up according to the staggered
57!>       grid. However, this implies no error since staggered velocity components
58!>       are zero at the walls and inside buildings.
59!------------------------------------------------------------------------------!
60 SUBROUTINE flow_statistics
61
62
63    USE arrays_3d,                                                             &
64        ONLY:  ddzu, ddzw, e, heatflux_output_conversion, hyp, km, kh,         &
65               momentumflux_output_conversion, nc, nr, p, prho, prr, pt, q,    &
66               qc, ql, qr, rho_air, rho_air_zw, rho_ocean, s,                  &
67               sa, u, ug, v, vg, vpt, w, w_subs, waterflux_output_conversion,  &
68               zw, d_exner
69
70    USE basic_constants_and_equations_mod,                                     &
71        ONLY:  g, lv_d_cp
72
73    USE bulk_cloud_model_mod,                                                  &
74        ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert
75
76    USE chem_modules,                                                          &
77        ONLY:  max_pr_cs
78
79    USE control_parameters,                                                    &
80        ONLY:   air_chemistry, average_count_pr, cloud_droplets, do_sum,       &
81                dt_3d, humidity, initializing_actions, land_surface,           &
82                large_scale_forcing, large_scale_subsidence, max_pr_user,      &
83                message_string, neutral, ocean_mode, passive_scalar,           &
84                simulated_time, simulated_time_at_begin,                       &
85                use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
86                ws_scheme_mom, ws_scheme_sca, salsa, max_pr_salsa
87
88    USE cpulog,                                                                &
89        ONLY:   cpu_log, log_point
90
91    USE grid_variables,                                                        &
92        ONLY:   ddx, ddy
93       
94    USE indices,                                                               &
95        ONLY:   ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums,      &
96                ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzt, topo_min_level,     &
97                wall_flags_static_0
98       
99    USE kinds
100   
101    USE land_surface_model_mod,                                                &
102        ONLY:   m_soil_h, nzb_soil, nzt_soil, t_soil_h
103
104    USE lsf_nudging_mod,                                                       &
105        ONLY:   td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q, time_vert
106
107    USE module_interface,                                                      &
108        ONLY:  module_interface_statistics
109
110    USE netcdf_interface,                                                      &
111        ONLY:  dots_rad, dots_soil, dots_max
112
113    USE pegrid
114
115    USE radiation_model_mod,                                                   &
116        ONLY:  radiation, radiation_scheme,                                    &
117               rad_lw_in, rad_lw_out, rad_lw_cs_hr, rad_lw_hr,                 &
118               rad_sw_in, rad_sw_out, rad_sw_cs_hr, rad_sw_hr
119
120    USE statistics
121
122    USE surface_mod,                                                           &
123        ONLY :  surf_def_h, surf_lsm_h, surf_usm_h
124
125
126    IMPLICIT NONE
127
128    INTEGER(iwp) ::  i                   !<
129    INTEGER(iwp) ::  j                   !<
130    INTEGER(iwp) ::  k                   !<
131    INTEGER(iwp) ::  ki                  !<
132    INTEGER(iwp) ::  k_surface_level     !<
133    INTEGER(iwp) ::  m                   !< loop variable over all horizontal wall elements
134    INTEGER(iwp) ::  l                   !< loop variable over surface facing -- up- or downward-facing
135    INTEGER(iwp) ::  nt                  !<
136!$  INTEGER(iwp) ::  omp_get_thread_num  !<
137    INTEGER(iwp) ::  sr                  !<
138    INTEGER(iwp) ::  tn                  !<
139
140    LOGICAL ::  first  !<
141   
142    REAL(wp) ::  dptdz_threshold  !<
143    REAL(wp) ::  fac              !<
144    REAL(wp) ::  flag             !<
145    REAL(wp) ::  height           !<
146    REAL(wp) ::  pts              !<
147    REAL(wp) ::  sums_l_etot      !<
148    REAL(wp) ::  ust              !<
149    REAL(wp) ::  ust2             !<
150    REAL(wp) ::  u2               !<
151    REAL(wp) ::  vst              !<
152    REAL(wp) ::  vst2             !<
153    REAL(wp) ::  v2               !<
154    REAL(wp) ::  w2               !<
155   
156    REAL(wp) ::  dptdz(nzb+1:nzt+1)    !<
157    REAL(wp) ::  sums_ll(nzb:nzt+1,2)  !<
158
159    CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
160
161
162!
163!-- To be on the safe side, check whether flow_statistics has already been
164!-- called once after the current time step
165    IF ( flow_statistics_called )  THEN
166
167       message_string = 'flow_statistics is called two times within one ' // &
168                        'timestep'
169       CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
170
171    ENDIF
172
173!
174!-- Compute statistics for each (sub-)region
175    DO  sr = 0, statistic_regions
176
177!
178!--    Initialize (local) summation array
179       sums_l = 0.0_wp
180#ifdef _OPENACC
181       !$ACC KERNELS PRESENT(sums_l)
182       sums_l = 0.0_wp
183       !$ACC END KERNELS
184#endif
185
186!
187!--    Store sums that have been computed in other subroutines in summation
188!--    array
189       sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
190!--    WARNING: next line still has to be adjusted for OpenMP
191       sums_l(:,21,0) = sums_wsts_bc_l(:,sr) *                                 &
192                        heatflux_output_conversion  ! heat flux from advec_s_bc
193       sums_l(nzb+9,pr_palm,0)  = sums_divold_l(sr)  ! old divergence from pres
194       sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
195
196!
197!--    When calcuating horizontally-averaged total (resolved- plus subgrid-
198!--    scale) vertical fluxes and velocity variances by using commonly-
199!--    applied Reynolds-based methods ( e.g. <w'pt'> = (w-<w>)*(pt-<pt>) )
200!--    in combination with the 5th order advection scheme, pronounced
201!--    artificial kinks could be observed in the vertical profiles near the
202!--    surface. Please note: these kinks were not related to the model truth,
203!--    i.e. these kinks are just related to an evaluation problem.   
204!--    In order avoid these kinks, vertical fluxes and horizontal as well
205!--    vertical velocity variances are calculated directly within the advection
206!--    routines, according to the numerical discretization, to evaluate the
207!--    statistical quantities as they will appear within the prognostic
208!--    equations.
209!--    Copy the turbulent quantities, evaluated in the advection routines to
210!--    the local array sums_l() for further computations.
211       IF ( ws_scheme_mom .AND. sr == 0 )  THEN
212
213!
214!--       According to the Neumann bc for the horizontal velocity components,
215!--       the corresponding fluxes has to satisfiy the same bc.
216          IF ( ocean_mode )  THEN
217             sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
218             sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)
219          ENDIF
220
221          DO  i = 0, threads_per_task-1
222!
223!--          Swap the turbulent quantities evaluated in advec_ws.
224             sums_l(:,13,i) = sums_wsus_ws_l(:,i)                              &
225                              * momentumflux_output_conversion ! w*u*
226             sums_l(:,15,i) = sums_wsvs_ws_l(:,i)                              &
227                              * momentumflux_output_conversion ! w*v*
228             sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2
229             sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2
230             sums_l(:,32,i) = sums_ws2_ws_l(:,i)        ! w*2
231             sums_l(:,34,i) = sums_l(:,34,i) + 0.5_wp *                        & 
232                              ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +      &
233                                sums_ws2_ws_l(:,i) )    ! e*
234          ENDDO
235
236       ENDIF
237
238       IF ( ws_scheme_sca .AND. sr == 0 )  THEN
239
240          DO  i = 0, threads_per_task-1
241             sums_l(:,17,i)                        = sums_wspts_ws_l(:,i)      &
242                                           * heatflux_output_conversion  ! w*pt*
243             IF ( ocean_mode     ) sums_l(:,66,i)  = sums_wssas_ws_l(:,i) ! w*sa*
244             IF ( humidity       ) sums_l(:,49,i)  = sums_wsqs_ws_l(:,i)       &
245                                           * waterflux_output_conversion  ! w*q*
246             IF ( passive_scalar ) sums_l(:,114,i) = sums_wsss_ws_l(:,i)  ! w*s*
247          ENDDO
248
249       ENDIF
250!
251!--    Horizontally averaged profiles of horizontal velocities and temperature.
252!--    They must have been computed before, because they are already required
253!--    for other horizontal averages.
254       tn = 0
255       !$OMP PARALLEL PRIVATE( i, j, k, tn, flag )
256       !$ tn = omp_get_thread_num()
257       !$OMP DO
258       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k, flag) &
259       !$ACC PRESENT(wall_flags_static_0, u, v, pt, rmask, sums_l)
260       DO  i = nxl, nxr
261          DO  j =  nys, nyn
262             DO  k = nzb, nzt+1
263                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j,i), 22 ) )
264                !$ACC ATOMIC
265                sums_l(k,1,tn)  = sums_l(k,1,tn)  + u(k,j,i)  * rmask(j,i,sr)  &
266                                                              * flag
267                !$ACC ATOMIC
268                sums_l(k,2,tn)  = sums_l(k,2,tn)  + v(k,j,i)  * rmask(j,i,sr)  &
269                                                              * flag
270                !$ACC ATOMIC
271                sums_l(k,4,tn)  = sums_l(k,4,tn)  + pt(k,j,i) * rmask(j,i,sr)  &
272                                                              * flag
273             ENDDO
274          ENDDO
275       ENDDO
276       !$ACC UPDATE HOST(sums_l(:,1,tn), sums_l(:,2,tn), sums_l(:,4,tn))
277
278!
279!--    Horizontally averaged profile of salinity
280       IF ( ocean_mode )  THEN
281          !$OMP DO
282          DO  i = nxl, nxr
283             DO  j =  nys, nyn
284                DO  k = nzb, nzt+1
285                   sums_l(k,23,tn)  = sums_l(k,23,tn) + sa(k,j,i)              &
286                                    * rmask(j,i,sr)                            &
287                                    * MERGE( 1.0_wp, 0.0_wp,                   &
288                                             BTEST( wall_flags_static_0(k,j,i), 22 ) )
289                ENDDO
290             ENDDO
291          ENDDO
292       ENDIF
293
294!
295!--    Horizontally averaged profiles of virtual potential temperature,
296!--    total water content, water vapor mixing ratio and liquid water potential
297!--    temperature
298       IF ( humidity )  THEN
299          !$OMP DO
300          DO  i = nxl, nxr
301             DO  j =  nys, nyn
302                DO  k = nzb, nzt+1
303                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j,i), 22 ) )
304                   sums_l(k,44,tn)  = sums_l(k,44,tn) +                        &
305                                      vpt(k,j,i) * rmask(j,i,sr) * flag
306                   sums_l(k,41,tn)  = sums_l(k,41,tn) +                        &
307                                      q(k,j,i) * rmask(j,i,sr)   * flag
308                ENDDO
309             ENDDO
310          ENDDO
311          IF ( bulk_cloud_model )  THEN
312             !$OMP DO
313             DO  i = nxl, nxr
314                DO  j =  nys, nyn
315                   DO  k = nzb, nzt+1
316                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j,i), 22 ) )
317                      sums_l(k,42,tn) = sums_l(k,42,tn) +                      &
318                                      ( q(k,j,i) - ql(k,j,i) ) * rmask(j,i,sr) &
319                                                               * flag
320                      sums_l(k,43,tn) = sums_l(k,43,tn) + (                    &
321                                      pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i) &
322                                                          ) * rmask(j,i,sr)    &
323                                                            * flag
324                   ENDDO
325                ENDDO
326             ENDDO
327          ENDIF
328       ENDIF
329
330!
331!--    Horizontally averaged profiles of passive scalar
332       IF ( passive_scalar )  THEN
333          !$OMP DO
334          DO  i = nxl, nxr
335             DO  j =  nys, nyn
336                DO  k = nzb, nzt+1
337                   sums_l(k,115,tn)  = sums_l(k,115,tn) + s(k,j,i)             &
338                                    * rmask(j,i,sr)                            &
339                                    * MERGE( 1.0_wp, 0.0_wp,                   &
340                                             BTEST( wall_flags_static_0(k,j,i), 22 ) )
341                ENDDO
342             ENDDO
343          ENDDO
344       ENDIF
345       !$OMP END PARALLEL
346!
347!--    Summation of thread sums
348       IF ( threads_per_task > 1 )  THEN
349          DO  i = 1, threads_per_task-1
350             sums_l(:,1,0) = sums_l(:,1,0) + sums_l(:,1,i)
351             sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i)
352             sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i)
353             IF ( ocean_mode )  THEN
354                sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)
355             ENDIF
356             IF ( humidity )  THEN
357                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
358                sums_l(:,44,0) = sums_l(:,44,0) + sums_l(:,44,i)
359                IF ( bulk_cloud_model )  THEN
360                   sums_l(:,42,0) = sums_l(:,42,0) + sums_l(:,42,i)
361                   sums_l(:,43,0) = sums_l(:,43,0) + sums_l(:,43,i)
362                ENDIF
363             ENDIF
364             IF ( passive_scalar )  THEN
365                sums_l(:,115,0) = sums_l(:,115,0) + sums_l(:,115,i)
366             ENDIF
367          ENDDO
368       ENDIF
369
370#if defined( __parallel )
371!
372!--    Compute total sum from local sums
373       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
374       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL,  &
375                           MPI_SUM, comm2d, ierr )
376       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
377       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL,  &
378                           MPI_SUM, comm2d, ierr )
379       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
380       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL,  &
381                           MPI_SUM, comm2d, ierr )
382       IF ( ocean_mode )  THEN
383          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
384          CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb,       &
385                              MPI_REAL, MPI_SUM, comm2d, ierr )
386       ENDIF
387       IF ( humidity ) THEN
388          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
389          CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb,       &
390                              MPI_REAL, MPI_SUM, comm2d, ierr )
391          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
392          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb,       &
393                              MPI_REAL, MPI_SUM, comm2d, ierr )
394          IF ( bulk_cloud_model ) THEN
395             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
396             CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb,    &
397                                 MPI_REAL, MPI_SUM, comm2d, ierr )
398             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
399             CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb,    &
400                                 MPI_REAL, MPI_SUM, comm2d, ierr )
401          ENDIF
402       ENDIF
403
404       IF ( passive_scalar )  THEN
405          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
406          CALL MPI_ALLREDUCE( sums_l(nzb,115,0), sums(nzb,115), nzt+2-nzb,       &
407                              MPI_REAL, MPI_SUM, comm2d, ierr )
408       ENDIF
409#else
410       sums(:,1) = sums_l(:,1,0)
411       sums(:,2) = sums_l(:,2,0)
412       sums(:,4) = sums_l(:,4,0)
413       IF ( ocean_mode )  sums(:,23) = sums_l(:,23,0)
414       IF ( humidity ) THEN
415          sums(:,44) = sums_l(:,44,0)
416          sums(:,41) = sums_l(:,41,0)
417          IF ( bulk_cloud_model ) THEN
418             sums(:,42) = sums_l(:,42,0)
419             sums(:,43) = sums_l(:,43,0)
420          ENDIF
421       ENDIF
422       IF ( passive_scalar )  sums(:,115) = sums_l(:,115,0)
423#endif
424
425!
426!--    Final values are obtained by division by the total number of grid points
427!--    used for summation. After that store profiles.
428       sums(:,1) = sums(:,1) / ngp_2dh(sr)
429       sums(:,2) = sums(:,2) / ngp_2dh(sr)
430       sums(:,4) = sums(:,4) / ngp_2dh_s_inner(:,sr)
431       hom(:,1,1,sr) = sums(:,1)             ! u
432       hom(:,1,2,sr) = sums(:,2)             ! v
433       hom(:,1,4,sr) = sums(:,4)             ! pt
434       !$ACC UPDATE DEVICE(hom(:,1,1,sr), hom(:,1,2,sr), hom(:,1,4,sr))
435
436
437!
438!--    Salinity
439       IF ( ocean_mode )  THEN
440          sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr)
441          hom(:,1,23,sr) = sums(:,23)             ! sa
442       ENDIF
443
444!
445!--    Humidity and cloud parameters
446       IF ( humidity ) THEN
447          sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr)
448          sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
449          hom(:,1,44,sr) = sums(:,44)             ! vpt
450          hom(:,1,41,sr) = sums(:,41)             ! qv (q)
451          IF ( bulk_cloud_model ) THEN
452             sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr)
453             sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr)
454             hom(:,1,42,sr) = sums(:,42)             ! qv
455             hom(:,1,43,sr) = sums(:,43)             ! pt
456          ENDIF
457       ENDIF
458
459!
460!--    Passive scalar
461       IF ( passive_scalar )  hom(:,1,115,sr) = sums(:,115) /                  &
462            ngp_2dh_s_inner(:,sr)                    ! s
463
464!
465!--    Horizontally averaged profiles of the remaining prognostic variables,
466!--    variances, the total and the perturbation energy (single values in last
467!--    column of sums_l) and some diagnostic quantities.
468!--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
469!--    ----  speaking the following k-loop would have to be split up and
470!--          rearranged according to the staggered grid.
471!--          However, this implies no error since staggered velocity components
472!--          are zero at the walls and inside buildings.
473       tn = 0
474       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll,                          &
475       !$OMP                   sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2,  &
476       !$OMP                   w2, flag, m, ki, l )
477       !$ tn = omp_get_thread_num()
478       !$OMP DO
479       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, m) &
480       !$ACC PRIVATE(sums_l_etot, flag) &
481       !$ACC PRESENT(wall_flags_static_0, rmask, momentumflux_output_conversion) &
482       !$ACC PRESENT(hom(:,1,4,sr)) &
483       !$ACC PRESENT(e, u, v, w, km, kh, p, pt) &
484       !$ACC PRESENT(surf_def_h(0), surf_lsm_h, surf_usm_h) &
485       !$ACC PRESENT(sums_l)
486       DO  i = nxl, nxr
487          DO  j =  nys, nyn
488             sums_l_etot = 0.0_wp
489             DO  k = nzb, nzt+1
490                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j,i), 22 ) )
491!
492!--             Prognostic and diagnostic variables
493                !$ACC ATOMIC
494                sums_l(k,3,tn)  = sums_l(k,3,tn)  + w(k,j,i)  * rmask(j,i,sr)  &
495                                                              * flag
496                !$ACC ATOMIC
497                sums_l(k,8,tn)  = sums_l(k,8,tn)  + e(k,j,i)  * rmask(j,i,sr)  &
498                                                              * flag
499                !$ACC ATOMIC
500                sums_l(k,9,tn)  = sums_l(k,9,tn)  + km(k,j,i) * rmask(j,i,sr)  &
501                                                              * flag
502                !$ACC ATOMIC
503                sums_l(k,10,tn) = sums_l(k,10,tn) + kh(k,j,i) * rmask(j,i,sr)  &
504                                                              * flag
505                !$ACC ATOMIC
506                sums_l(k,40,tn) = sums_l(k,40,tn) + ( p(k,j,i)                 &
507                                         / momentumflux_output_conversion(k) ) &
508                                                              * flag
509
510                !$ACC ATOMIC
511                sums_l(k,33,tn) = sums_l(k,33,tn) + &
512                                  ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr)&
513                                                                 * flag
514#ifndef _OPENACC
515                IF ( humidity )  THEN
516                   sums_l(k,70,tn) = sums_l(k,70,tn) + &
517                                  ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr)&
518                                                                 * flag
519                ENDIF
520                IF ( passive_scalar )  THEN
521                   sums_l(k,116,tn) = sums_l(k,116,tn) + &
522                                  ( s(k,j,i)-hom(k,1,115,sr) )**2 * rmask(j,i,sr)&
523                                                                  * flag
524                ENDIF
525#endif
526!
527!--             Higher moments
528!--             (Computation of the skewness of w further below)
529                !$ACC ATOMIC
530                sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i)**3 * rmask(j,i,sr) &
531                                                                * flag
532
533                sums_l_etot  = sums_l_etot + &
534                                        0.5_wp * ( u(k,j,i)**2 + v(k,j,i)**2 +  &
535                                        w(k,j,i)**2 )            * rmask(j,i,sr)&
536                                                                 * flag
537             ENDDO
538!
539!--          Total and perturbation energy for the total domain (being
540!--          collected in the last column of sums_l). Summation of these
541!--          quantities is seperated from the previous loop in order to
542!--          allow vectorization of that loop.
543             !$ACC ATOMIC
544             sums_l(nzb+4,pr_palm,tn) = sums_l(nzb+4,pr_palm,tn) + sums_l_etot
545!
546!--          2D-arrays (being collected in the last column of sums_l)
547             IF ( surf_def_h(0)%end_index(j,i) >=                              &
548                  surf_def_h(0)%start_index(j,i) )  THEN
549                m = surf_def_h(0)%start_index(j,i)
550                !$ACC ATOMIC
551                sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +            &
552                                        surf_def_h(0)%us(m)   * rmask(j,i,sr)
553                !$ACC ATOMIC
554                sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +          &
555                                        surf_def_h(0)%usws(m) * rmask(j,i,sr)
556                !$ACC ATOMIC
557                sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +          &
558                                        surf_def_h(0)%vsws(m) * rmask(j,i,sr)
559                !$ACC ATOMIC
560                sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +          &
561                                        surf_def_h(0)%ts(m)   * rmask(j,i,sr)
562#ifndef _OPENACC
563                IF ( humidity )  THEN
564                   sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +     &
565                                            surf_def_h(0)%qs(m)   * rmask(j,i,sr)
566                ENDIF
567                IF ( passive_scalar )  THEN
568                   sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) +     &
569                                            surf_def_h(0)%ss(m)   * rmask(j,i,sr)
570                ENDIF
571#endif
572!
573!--             Summation of surface temperature.
574                !$ACC ATOMIC
575                sums_l(nzb+14,pr_palm,tn) = sums_l(nzb+14,pr_palm,tn)   +      &
576                                            surf_def_h(0)%pt_surface(m) *      &
577                                            rmask(j,i,sr)
578             ENDIF
579             IF ( surf_lsm_h%end_index(j,i) >= surf_lsm_h%start_index(j,i) )  THEN
580                m = surf_lsm_h%start_index(j,i)
581                !$ACC ATOMIC
582                sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +            &
583                                        surf_lsm_h%us(m)   * rmask(j,i,sr)
584                !$ACC ATOMIC
585                sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +          &
586                                        surf_lsm_h%usws(m) * rmask(j,i,sr)
587                !$ACC ATOMIC
588                sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +          &
589                                        surf_lsm_h%vsws(m) * rmask(j,i,sr)
590                !$ACC ATOMIC
591                sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +          &
592                                        surf_lsm_h%ts(m)   * rmask(j,i,sr)
593#ifndef _OPENACC
594                IF ( humidity )  THEN
595                   sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +     &
596                                            surf_lsm_h%qs(m)   * rmask(j,i,sr)
597                ENDIF
598                IF ( passive_scalar )  THEN
599                   sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) +     &
600                                            surf_lsm_h%ss(m)   * rmask(j,i,sr)
601                ENDIF
602#endif
603!
604!--             Summation of surface temperature.
605                !$ACC ATOMIC
606                sums_l(nzb+14,pr_palm,tn) = sums_l(nzb+14,pr_palm,tn)   +      &
607                                            surf_lsm_h%pt_surface(m)    *      &
608                                            rmask(j,i,sr)
609             ENDIF
610             IF ( surf_usm_h%end_index(j,i) >= surf_usm_h%start_index(j,i) )  THEN
611                m = surf_usm_h%start_index(j,i)
612                !$ACC ATOMIC
613                sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +            &
614                                        surf_usm_h%us(m)   * rmask(j,i,sr)
615                !$ACC ATOMIC
616                sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +          &
617                                        surf_usm_h%usws(m) * rmask(j,i,sr)
618                !$ACC ATOMIC
619                sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +          &
620                                        surf_usm_h%vsws(m) * rmask(j,i,sr)
621                !$ACC ATOMIC
622                sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +          &
623                                        surf_usm_h%ts(m)   * rmask(j,i,sr)
624#ifndef _OPENACC
625                IF ( humidity )  THEN
626                   sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +     &
627                                            surf_usm_h%qs(m)   * rmask(j,i,sr)
628                ENDIF
629                IF ( passive_scalar )  THEN
630                   sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) +     &
631                                            surf_usm_h%ss(m)   * rmask(j,i,sr)
632                ENDIF
633#endif
634!
635!--             Summation of surface temperature.
636                !$ACC ATOMIC
637                sums_l(nzb+14,pr_palm,tn) = sums_l(nzb+14,pr_palm,tn)   +      &
638                                            surf_usm_h%pt_surface(m)    *      &
639                                            rmask(j,i,sr)
640             ENDIF
641          ENDDO
642       ENDDO
643       !$ACC UPDATE &
644       !$ACC HOST(sums_l(:,3,tn), sums_l(:,8,tn), sums_l(:,9,tn)) &
645       !$ACC HOST(sums_l(:,10,tn), sums_l(:,40,tn), sums_l(:,33,tn)) &
646       !$ACC HOST(sums_l(:,38,tn)) &
647       !$ACC HOST(sums_l(nzb:nzb+4,pr_palm,tn), sums_l(nzb+14:nzb+14,pr_palm,tn))
648
649!
650!--    Computation of statistics when ws-scheme is not used. Else these
651!--    quantities are evaluated in the advection routines.
652       IF ( .NOT. ws_scheme_mom .OR. sr /= 0 .OR. simulated_time == 0.0_wp )   &
653       THEN
654          !$OMP DO
655          DO  i = nxl, nxr
656             DO  j =  nys, nyn
657                DO  k = nzb, nzt+1
658                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j,i), 22 ) )
659
660                   u2   = u(k,j,i)**2
661                   v2   = v(k,j,i)**2
662                   w2   = w(k,j,i)**2
663                   ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
664                   vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
665
666                   sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr)    &
667                                                            * flag
668                   sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr)    &
669                                                            * flag
670                   sums_l(k,32,tn) = sums_l(k,32,tn) + w2   * rmask(j,i,sr)    &
671                                                            * flag
672!
673!--                Perturbation energy
674
675                   sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5_wp *                &
676                                  ( ust2 + vst2 + w2 )      * rmask(j,i,sr)    &
677                                                            * flag
678                ENDDO
679             ENDDO
680          ENDDO
681       ENDIF
682!
683!--    Computaion of domain-averaged perturbation energy. Please note,
684!--    to prevent that perturbation energy is larger (even if only slightly)
685!--    than the total kinetic energy, calculation is based on deviations from
686!--    the horizontal mean, instead of spatial descretization of the advection
687!--    term.
688       !$OMP DO
689       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k, flag, w2, ust2, vst2) &
690       !$ACC PRESENT(wall_flags_static_0, u, v, w, rmask, hom(:,1,1:2,sr)) &
691       !$ACC PRESENT(sums_l)
692       DO  i = nxl, nxr
693          DO  j =  nys, nyn
694             DO  k = nzb, nzt+1
695                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j,i), 22 ) )
696
697                w2   = w(k,j,i)**2
698                ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
699                vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
700                w2   = w(k,j,i)**2
701
702                !$ACC ATOMIC
703                sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn)            &
704                                 + 0.5_wp * ( ust2 + vst2 + w2 )               &
705                                 * rmask(j,i,sr)                               &
706                                 * flag
707             ENDDO
708          ENDDO
709       ENDDO
710       !$ACC UPDATE HOST(sums_l(nzb+5:nzb+5,pr_palm,tn))
711
712!
713!--    Horizontally averaged profiles of the vertical fluxes
714
715       !$OMP DO
716       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, l, m) &
717       !$ACC PRIVATE(ki, flag, ust, vst, pts) &
718       !$ACC PRESENT(kh, km, u, v, w, pt) &
719       !$ACC PRESENT(wall_flags_static_0, rmask, ddzu, rho_air_zw, hom(:,1,1:4,sr)) &
720       !$ACC PRESENT(heatflux_output_conversion, momentumflux_output_conversion) &
721       !$ACC PRESENT(surf_def_h(0:2), surf_lsm_h, surf_usm_h) &
722       !$ACC PRESENT(sums_l)
723       DO  i = nxl, nxr
724          DO  j = nys, nyn
725!
726!--          Subgridscale fluxes (without Prandtl layer from k=nzb,
727!--          oterwise from k=nzb+1)
728!--          NOTE: for simplicity, nzb_diff_s_inner is used below, although
729!--          ----  strictly speaking the following k-loop would have to be
730!--                split up according to the staggered grid.
731!--                However, this implies no error since staggered velocity
732!--                components are zero at the walls and inside buildings.
733!--          Flag 23 is used to mask surface fluxes as well as model-top fluxes,
734!--          which are added further below.
735             DO  k = nzb, nzt
736                flag = MERGE( 1.0_wp, 0.0_wp,                                  &
737                              BTEST( wall_flags_static_0(k,j,i), 23 ) ) *             &
738                       MERGE( 1.0_wp, 0.0_wp,                                  &
739                              BTEST( wall_flags_static_0(k,j,i), 9  ) )
740!
741!--             Momentum flux w"u"
742                !$ACC ATOMIC
743                sums_l(k,12,tn) = sums_l(k,12,tn) - 0.25_wp * (                &
744                               km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
745                                                           ) * (               &
746                                   ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
747                                 + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
748                                                           ) * rmask(j,i,sr)   &
749                                         * rho_air_zw(k)                       &
750                                         * momentumflux_output_conversion(k)   &
751                                         * flag
752!
753!--             Momentum flux w"v"
754                !$ACC ATOMIC
755                sums_l(k,14,tn) = sums_l(k,14,tn) - 0.25_wp * (                &
756                               km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
757                                                           ) * (               &
758                                   ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
759                                 + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
760                                                           ) * rmask(j,i,sr)   &
761                                         * rho_air_zw(k)                       &
762                                         * momentumflux_output_conversion(k)   &
763                                         * flag
764!
765!--             Heat flux w"pt"
766                !$ACC ATOMIC
767                sums_l(k,16,tn) = sums_l(k,16,tn)                              &
768                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
769                                               * ( pt(k+1,j,i) - pt(k,j,i) )   &
770                                               * rho_air_zw(k)                 &
771                                               * heatflux_output_conversion(k) &
772                                               * ddzu(k+1) * rmask(j,i,sr)     &
773                                               * flag
774
775!
776!--             Salinity flux w"sa"
777#ifndef _OPENACC
778                IF ( ocean_mode )  THEN
779                   sums_l(k,65,tn) = sums_l(k,65,tn)                           &
780                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
781                                               * ( sa(k+1,j,i) - sa(k,j,i) )   &
782                                               * ddzu(k+1) * rmask(j,i,sr)     &
783                                               * flag
784                ENDIF
785
786!
787!--             Buoyancy flux, water flux (humidity flux) w"q"
788                IF ( humidity ) THEN
789                   sums_l(k,45,tn) = sums_l(k,45,tn)                           &
790                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
791                                               * ( vpt(k+1,j,i) - vpt(k,j,i) ) &
792                                               * rho_air_zw(k)                 &
793                                               * heatflux_output_conversion(k) &
794                                               * ddzu(k+1) * rmask(j,i,sr) * flag
795                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
796                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
797                                               * ( q(k+1,j,i) - q(k,j,i) )     &
798                                               * rho_air_zw(k)                 &
799                                               * waterflux_output_conversion(k)&
800                                               * ddzu(k+1) * rmask(j,i,sr) * flag
801
802                   IF ( bulk_cloud_model ) THEN
803                      sums_l(k,51,tn) = sums_l(k,51,tn)                        &
804                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
805                                               * ( ( q(k+1,j,i) - ql(k+1,j,i) )&
806                                                - ( q(k,j,i) - ql(k,j,i) ) )   &
807                                               * rho_air_zw(k)                 &
808                                               * waterflux_output_conversion(k)&
809                                               * ddzu(k+1) * rmask(j,i,sr) * flag
810                   ENDIF
811                ENDIF
812
813!
814!--             Passive scalar flux
815                IF ( passive_scalar )  THEN
816                   sums_l(k,117,tn) = sums_l(k,117,tn)                         &
817                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
818                                                  * ( s(k+1,j,i) - s(k,j,i) )  &
819                                                  * ddzu(k+1) * rmask(j,i,sr)  &
820                                                              * flag
821                ENDIF
822#endif
823
824             ENDDO
825
826!
827!--          Subgridscale fluxes in the Prandtl layer
828             IF ( use_surface_fluxes )  THEN
829                DO  l = 0, 1
830                   ! The original code using MERGE doesn't work with the PGI
831                   ! compiler when running on the GPU.
832                   ! This is submitted as a compiler Bug in PGI ticket TPR#26718
833                   ! ki = MERGE( -1, 0, l == 0 )
834                   ki = -1 + l
835                   IF ( surf_def_h(l)%ns >= 1 )  THEN
836                      DO  m = surf_def_h(l)%start_index(j,i),                  &
837                              surf_def_h(l)%end_index(j,i)
838                         k = surf_def_h(l)%k(m)
839
840                         !$ACC ATOMIC
841                         sums_l(k+ki,12,tn) = sums_l(k+ki,12,tn) + &
842                                    momentumflux_output_conversion(k+ki) * &
843                                    surf_def_h(l)%usws(m) * rmask(j,i,sr)     ! w"u"
844                         !$ACC ATOMIC
845                         sums_l(k+ki,14,tn) = sums_l(k+ki,14,tn) + &
846                                    momentumflux_output_conversion(k+ki) * &
847                                    surf_def_h(l)%vsws(m) * rmask(j,i,sr)     ! w"v"
848                         !$ACC ATOMIC
849                         sums_l(k+ki,16,tn) = sums_l(k+ki,16,tn) + &
850                                    heatflux_output_conversion(k+ki) * &
851                                    surf_def_h(l)%shf(m)  * rmask(j,i,sr)     ! w"pt"
852#if 0
853                         sums_l(k+ki,58,tn) = sums_l(k+ki,58,tn) + &
854                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
855                         sums_l(k+ki,61,tn) = sums_l(k+ki,61,tn) + &
856                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
857#endif
858#ifndef _OPENACC
859                         IF ( ocean_mode )  THEN
860                            sums_l(k+ki,65,tn) = sums_l(k+ki,65,tn) + &
861                                       surf_def_h(l)%sasws(m) * rmask(j,i,sr)  ! w"sa"
862                         ENDIF
863                         IF ( humidity )  THEN
864                            sums_l(k+ki,48,tn) = sums_l(k+ki,48,tn) +                     &
865                                       waterflux_output_conversion(k+ki) *      &
866                                       surf_def_h(l)%qsws(m) * rmask(j,i,sr)  ! w"q" (w"qv")
867                            sums_l(k+ki,45,tn) = sums_l(k+ki,45,tn) + (                   &
868                                       ( 1.0_wp + 0.61_wp * q(k+ki,j,i) ) *     &
869                                       surf_def_h(l)%shf(m) + 0.61_wp * pt(k+ki,j,i) *      &
870                                                  surf_def_h(l)%qsws(m) )                  &
871                                       * heatflux_output_conversion(k+ki)
872                            IF ( cloud_droplets )  THEN
873                               sums_l(k+ki,45,tn) = sums_l(k+ki,45,tn) + (                &
874                                         ( 1.0_wp + 0.61_wp * q(k+ki,j,i) -     &
875                                           ql(k+ki,j,i) ) * surf_def_h(l)%shf(m) +          &
876                                           0.61_wp * pt(k+ki,j,i) * surf_def_h(l)%qsws(m) ) &
877                                          * heatflux_output_conversion(k+ki)
878                            ENDIF
879                            IF ( bulk_cloud_model )  THEN
880!
881!--                            Formula does not work if ql(k+ki) /= 0.0
882                               sums_l(k+ki,51,tn) = sums_l(k+ki,51,tn) +                  &
883                                          waterflux_output_conversion(k+ki) *   &
884                                          surf_def_h(l)%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
885                            ENDIF
886                         ENDIF
887                         IF ( passive_scalar )  THEN
888                            sums_l(k+ki,117,tn) = sums_l(k+ki,117,tn) +                     &
889                                        surf_def_h(l)%ssws(m) * rmask(j,i,sr) ! w"s"
890                         ENDIF
891#endif
892
893                      ENDDO
894
895                   ENDIF
896                ENDDO
897                IF ( surf_lsm_h%end_index(j,i) >=                              &
898                     surf_lsm_h%start_index(j,i) )  THEN
899                   m = surf_lsm_h%start_index(j,i)
900                   !$ACC ATOMIC
901                   sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
902                                    momentumflux_output_conversion(nzb) * &
903                                    surf_lsm_h%usws(m) * rmask(j,i,sr)     ! w"u"
904                   !$ACC ATOMIC
905                   sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
906                                    momentumflux_output_conversion(nzb) * &
907                                    surf_lsm_h%vsws(m) * rmask(j,i,sr)     ! w"v"
908                   !$ACC ATOMIC
909                   sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
910                                    heatflux_output_conversion(nzb) * &
911                                    surf_lsm_h%shf(m)  * rmask(j,i,sr)     ! w"pt"
912#if 0
913                   sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
914                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
915                   sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
916                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
917#endif
918#ifndef _OPENACC
919                   IF ( ocean_mode )  THEN
920                      sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
921                                       surf_lsm_h%sasws(m) * rmask(j,i,sr)  ! w"sa"
922                   ENDIF
923                   IF ( humidity )  THEN
924                      sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
925                                       waterflux_output_conversion(nzb) *      &
926                                       surf_lsm_h%qsws(m) * rmask(j,i,sr)  ! w"q" (w"qv")
927                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                   &
928                                       ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) *     &
929                                       surf_lsm_h%shf(m) + 0.61_wp * pt(nzb,j,i) *      &
930                                                  surf_lsm_h%qsws(m) )                  &
931                                       * heatflux_output_conversion(nzb)
932                      IF ( cloud_droplets )  THEN
933                         sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                &
934                                         ( 1.0_wp + 0.61_wp * q(nzb,j,i) -     &
935                                           ql(nzb,j,i) ) * surf_lsm_h%shf(m) +          &
936                                           0.61_wp * pt(nzb,j,i) * surf_lsm_h%qsws(m) ) &
937                                          * heatflux_output_conversion(nzb)
938                      ENDIF
939                      IF ( bulk_cloud_model )  THEN
940!
941!--                      Formula does not work if ql(nzb) /= 0.0
942                         sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                  &
943                                          waterflux_output_conversion(nzb) *   &
944                                          surf_lsm_h%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
945                      ENDIF
946                   ENDIF
947                   IF ( passive_scalar )  THEN
948                      sums_l(nzb,117,tn) = sums_l(nzb,117,tn) +                     &
949                                        surf_lsm_h%ssws(m) * rmask(j,i,sr) ! w"s"
950                   ENDIF
951#endif
952
953                ENDIF
954                IF ( surf_usm_h%end_index(j,i) >=                              &
955                     surf_usm_h%start_index(j,i) )  THEN
956                   m = surf_usm_h%start_index(j,i)
957                   !$ACC ATOMIC
958                   sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
959                                    momentumflux_output_conversion(nzb) * &
960                                    surf_usm_h%usws(m) * rmask(j,i,sr)     ! w"u"
961                   !$ACC ATOMIC
962                   sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
963                                    momentumflux_output_conversion(nzb) * &
964                                    surf_usm_h%vsws(m) * rmask(j,i,sr)     ! w"v"
965                   !$ACC ATOMIC
966                   sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
967                                    heatflux_output_conversion(nzb) * &
968                                    surf_usm_h%shf(m)  * rmask(j,i,sr)     ! w"pt"
969#if 0
970                   sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
971                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
972                   sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
973                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
974#endif
975#ifndef _OPENACC
976                   IF ( ocean_mode )  THEN
977                      sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
978                                       surf_usm_h%sasws(m) * rmask(j,i,sr)  ! w"sa"
979                   ENDIF
980                   IF ( humidity )  THEN
981                      sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
982                                       waterflux_output_conversion(nzb) *      &
983                                       surf_usm_h%qsws(m) * rmask(j,i,sr)  ! w"q" (w"qv")
984                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                   &
985                                       ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) *     &
986                                       surf_usm_h%shf(m) + 0.61_wp * pt(nzb,j,i) *      &
987                                                  surf_usm_h%qsws(m) )                  &
988                                       * heatflux_output_conversion(nzb)
989                      IF ( cloud_droplets )  THEN
990                         sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                &
991                                         ( 1.0_wp + 0.61_wp * q(nzb,j,i) -     &
992                                           ql(nzb,j,i) ) * surf_usm_h%shf(m) +          &
993                                           0.61_wp * pt(nzb,j,i) * surf_usm_h%qsws(m) ) &
994                                          * heatflux_output_conversion(nzb)
995                      ENDIF
996                      IF ( bulk_cloud_model )  THEN
997!
998!--                      Formula does not work if ql(nzb) /= 0.0
999                         sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                  &
1000                                          waterflux_output_conversion(nzb) *   &
1001                                          surf_usm_h%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
1002                      ENDIF
1003                   ENDIF
1004                   IF ( passive_scalar )  THEN
1005                      sums_l(nzb,117,tn) = sums_l(nzb,117,tn) +                     &
1006                                        surf_usm_h%ssws(m) * rmask(j,i,sr) ! w"s"
1007                   ENDIF
1008#endif
1009
1010                ENDIF
1011
1012             ENDIF
1013
1014#ifndef _OPENACC
1015             IF ( .NOT. neutral )  THEN
1016                IF ( surf_def_h(0)%end_index(j,i) >=                           &
1017                     surf_def_h(0)%start_index(j,i) )  THEN
1018                   m = surf_def_h(0)%start_index(j,i)
1019                   sums_l(nzb,112,tn) = sums_l(nzb,112,tn) +                   &
1020                                        surf_def_h(0)%ol(m)  * rmask(j,i,sr) ! L
1021                ENDIF
1022                IF ( surf_lsm_h%end_index(j,i) >=                              &
1023                     surf_lsm_h%start_index(j,i) )  THEN
1024                   m = surf_lsm_h%start_index(j,i)
1025                   sums_l(nzb,112,tn) = sums_l(nzb,112,tn) +                   &
1026                                        surf_lsm_h%ol(m)  * rmask(j,i,sr) ! L
1027                ENDIF
1028                IF ( surf_usm_h%end_index(j,i) >=                              &
1029                     surf_usm_h%start_index(j,i) )  THEN
1030                   m = surf_usm_h%start_index(j,i)
1031                   sums_l(nzb,112,tn) = sums_l(nzb,112,tn) +                   &
1032                                        surf_usm_h%ol(m)  * rmask(j,i,sr) ! L
1033                ENDIF
1034             ENDIF
1035
1036             IF ( radiation )  THEN
1037                IF ( surf_def_h(0)%end_index(j,i) >=                           &
1038                     surf_def_h(0)%start_index(j,i) )  THEN
1039                   m = surf_def_h(0)%start_index(j,i)
1040                   sums_l(nzb,99,tn)  = sums_l(nzb,99,tn)  +                   &
1041                                        surf_def_h(0)%rad_net(m) * rmask(j,i,sr)
1042                   sums_l(nzb,100,tn) = sums_l(nzb,100,tn)  +                  &
1043                                        surf_def_h(0)%rad_lw_in(m) * rmask(j,i,sr)
1044                   sums_l(nzb,101,tn) = sums_l(nzb,101,tn)  +                  &
1045                                        surf_def_h(0)%rad_lw_out(m) * rmask(j,i,sr)
1046                   sums_l(nzb,102,tn) = sums_l(nzb,102,tn)  +                  &
1047                                        surf_def_h(0)%rad_sw_in(m) * rmask(j,i,sr)
1048                   sums_l(nzb,103,tn) = sums_l(nzb,103,tn)  +                  &
1049                                        surf_def_h(0)%rad_sw_out(m) * rmask(j,i,sr)
1050                ENDIF
1051                IF ( surf_lsm_h%end_index(j,i) >=                              &
1052                     surf_lsm_h%start_index(j,i) )  THEN
1053                   m = surf_lsm_h%start_index(j,i)
1054                   sums_l(nzb,99,tn)  = sums_l(nzb,99,tn)  +                   &
1055                                        surf_lsm_h%rad_net(m) * rmask(j,i,sr)
1056                   sums_l(nzb,100,tn) = sums_l(nzb,100,tn)  +                  &
1057                                        surf_lsm_h%rad_lw_in(m) * rmask(j,i,sr)
1058                   sums_l(nzb,101,tn) = sums_l(nzb,101,tn)  +                  &
1059                                        surf_lsm_h%rad_lw_out(m) * rmask(j,i,sr)
1060                   sums_l(nzb,102,tn) = sums_l(nzb,102,tn)  +                  &
1061                                        surf_lsm_h%rad_sw_in(m) * rmask(j,i,sr)
1062                   sums_l(nzb,103,tn) = sums_l(nzb,103,tn)  +                  &
1063                                        surf_lsm_h%rad_sw_out(m) * rmask(j,i,sr)
1064                ENDIF
1065                IF ( surf_usm_h%end_index(j,i) >=                              &
1066                     surf_usm_h%start_index(j,i) )  THEN
1067                   m = surf_usm_h%start_index(j,i)
1068                   sums_l(nzb,99,tn)  = sums_l(nzb,99,tn)  +                   &
1069                                        surf_usm_h%rad_net(m) * rmask(j,i,sr)
1070                   sums_l(nzb,100,tn) = sums_l(nzb,100,tn)  +                  &
1071                                        surf_usm_h%rad_lw_in(m) * rmask(j,i,sr)
1072                   sums_l(nzb,101,tn) = sums_l(nzb,101,tn)  +                  &
1073                                        surf_usm_h%rad_lw_out(m) * rmask(j,i,sr)
1074                   sums_l(nzb,102,tn) = sums_l(nzb,102,tn)  +                  &
1075                                        surf_usm_h%rad_sw_in(m) * rmask(j,i,sr)
1076                   sums_l(nzb,103,tn) = sums_l(nzb,103,tn)  +                  &
1077                                        surf_usm_h%rad_sw_out(m) * rmask(j,i,sr)
1078                ENDIF
1079
1080#if defined ( __rrtmg )
1081                IF ( radiation_scheme == 'rrtmg' )  THEN
1082
1083                   IF ( surf_def_h(0)%end_index(j,i) >=                        &
1084                        surf_def_h(0)%start_index(j,i) )  THEN
1085                      m = surf_def_h(0)%start_index(j,i)
1086                      sums_l(nzb,108,tn)  = sums_l(nzb,108,tn)  +              &
1087                                   surf_def_h(0)%rrtm_aldif(0,m) * rmask(j,i,sr)
1088                      sums_l(nzb,109,tn) = sums_l(nzb,109,tn)  +               &
1089                                   surf_def_h(0)%rrtm_aldir(0,m) * rmask(j,i,sr)
1090                      sums_l(nzb,110,tn) = sums_l(nzb,110,tn)  +               &
1091                                   surf_def_h(0)%rrtm_asdif(0,m) * rmask(j,i,sr)
1092                      sums_l(nzb,111,tn) = sums_l(nzb,111,tn)  +               &
1093                                   surf_def_h(0)%rrtm_asdir(0,m) * rmask(j,i,sr)
1094                   ENDIF
1095                   IF ( surf_lsm_h%end_index(j,i) >=                           &
1096                        surf_lsm_h%start_index(j,i) )  THEN
1097                      m = surf_lsm_h%start_index(j,i)
1098                      sums_l(nzb,108,tn)  = sums_l(nzb,108,tn)  +              &
1099                               SUM( surf_lsm_h%frac(:,m) *                     &
1100                                    surf_lsm_h%rrtm_aldif(:,m) ) * rmask(j,i,sr)
1101                      sums_l(nzb,109,tn) = sums_l(nzb,109,tn)  +               &
1102                               SUM( surf_lsm_h%frac(:,m) *                     &
1103                                    surf_lsm_h%rrtm_aldir(:,m) ) * rmask(j,i,sr)
1104                      sums_l(nzb,110,tn) = sums_l(nzb,110,tn)  +               &
1105                               SUM( surf_lsm_h%frac(:,m) *                     &
1106                                    surf_lsm_h%rrtm_asdif(:,m) ) * rmask(j,i,sr)
1107                      sums_l(nzb,111,tn) = sums_l(nzb,111,tn)  +               &
1108                               SUM( surf_lsm_h%frac(:,m) *                     &
1109                                    surf_lsm_h%rrtm_asdir(:,m) ) * rmask(j,i,sr)
1110                   ENDIF
1111                   IF ( surf_usm_h%end_index(j,i) >=                           &
1112                        surf_usm_h%start_index(j,i) )  THEN
1113                      m = surf_usm_h%start_index(j,i)
1114                      sums_l(nzb,108,tn)  = sums_l(nzb,108,tn)  +              &
1115                               SUM( surf_usm_h%frac(:,m) *                     &
1116                                    surf_usm_h%rrtm_aldif(:,m) ) * rmask(j,i,sr)
1117                      sums_l(nzb,109,tn) = sums_l(nzb,109,tn)  +               &
1118                               SUM( surf_usm_h%frac(:,m) *                     &
1119                                    surf_usm_h%rrtm_aldir(:,m) ) * rmask(j,i,sr)
1120                      sums_l(nzb,110,tn) = sums_l(nzb,110,tn)  +               &
1121                               SUM( surf_usm_h%frac(:,m) *                     &
1122                                    surf_usm_h%rrtm_asdif(:,m) ) * rmask(j,i,sr)
1123                      sums_l(nzb,111,tn) = sums_l(nzb,111,tn)  +               &
1124                               SUM( surf_usm_h%frac(:,m) *                     &
1125                                    surf_usm_h%rrtm_asdir(:,m) ) * rmask(j,i,sr)
1126                   ENDIF
1127
1128                ENDIF
1129#endif
1130             ENDIF
1131#endif
1132!
1133!--          Subgridscale fluxes at the top surface
1134             IF ( use_top_fluxes )  THEN
1135                m = surf_def_h(2)%start_index(j,i)
1136                !$ACC ATOMIC
1137                sums_l(nzt,12,tn) = sums_l(nzt,12,tn) + &
1138                                    momentumflux_output_conversion(nzt) * &
1139                                    surf_def_h(2)%usws(m) * rmask(j,i,sr)    ! w"u"
1140                !$ACC ATOMIC
1141                sums_l(nzt+1,12,tn) = sums_l(nzt+1,12,tn) + &
1142                                    momentumflux_output_conversion(nzt+1) * &
1143                                    surf_def_h(2)%usws(m) * rmask(j,i,sr)    ! w"u"
1144                !$ACC ATOMIC
1145                sums_l(nzt,14,tn) = sums_l(nzt,14,tn) + &
1146                                    momentumflux_output_conversion(nzt) * &
1147                                    surf_def_h(2)%vsws(m) * rmask(j,i,sr)    ! w"v"
1148                !$ACC ATOMIC
1149                sums_l(nzt+1,14,tn) = sums_l(nzt+1,14,tn) + &
1150                                    momentumflux_output_conversion(nzt+1) * &
1151                                    surf_def_h(2)%vsws(m) * rmask(j,i,sr)    ! w"v"
1152                !$ACC ATOMIC
1153                sums_l(nzt,16,tn) = sums_l(nzt,16,tn) + &
1154                                    heatflux_output_conversion(nzt) * &
1155                                    surf_def_h(2)%shf(m)  * rmask(j,i,sr)   ! w"pt"
1156                !$ACC ATOMIC
1157                sums_l(nzt+1,16,tn) = sums_l(nzt+1,16,tn) + &
1158                                    heatflux_output_conversion(nzt+1) * &
1159                                    surf_def_h(2)%shf(m)  * rmask(j,i,sr)   ! w"pt"
1160#if 0
1161                sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + &
1162                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
1163                sums_l(nzt:nzt+1,61,tn) = sums_l(nzt:nzt+1,61,tn) + &
1164                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
1165#endif
1166#ifndef _OPENACC
1167                IF ( ocean_mode )  THEN
1168                   sums_l(nzt,65,tn) = sums_l(nzt,65,tn) + &
1169                                       surf_def_h(2)%sasws(m) * rmask(j,i,sr)  ! w"sa"
1170                ENDIF
1171                IF ( humidity )  THEN
1172                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) +                     &
1173                                       waterflux_output_conversion(nzt) *      &
1174                                       surf_def_h(2)%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
1175                   sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                   &
1176                                       ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) *     &
1177                                       surf_def_h(2)%shf(m) +                  &
1178                                       0.61_wp * pt(nzt,j,i) *    &
1179                                       surf_def_h(2)%qsws(m) )      &
1180                                       * heatflux_output_conversion(nzt)
1181                   IF ( cloud_droplets )  THEN
1182                      sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                &
1183                                          ( 1.0_wp + 0.61_wp * q(nzt,j,i) -    &
1184                                            ql(nzt,j,i) ) *                    &
1185                                            surf_def_h(2)%shf(m) +             &
1186                                           0.61_wp * pt(nzt,j,i) *             &
1187                                           surf_def_h(2)%qsws(m) )&
1188                                           * heatflux_output_conversion(nzt)
1189                   ENDIF
1190                   IF ( bulk_cloud_model )  THEN
1191!
1192!--                   Formula does not work if ql(nzb) /= 0.0
1193                      sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + &   ! w"q" (w"qv")
1194                                          waterflux_output_conversion(nzt) *   &
1195                                          surf_def_h(2)%qsws(m) * rmask(j,i,sr)
1196                   ENDIF
1197                ENDIF
1198                IF ( passive_scalar )  THEN
1199                   sums_l(nzt,117,tn) = sums_l(nzt,117,tn) + &
1200                                        surf_def_h(2)%ssws(m) * rmask(j,i,sr) ! w"s"
1201                ENDIF
1202#endif
1203             ENDIF
1204
1205!
1206!--          Resolved fluxes (can be computed for all horizontal points)
1207!--          NOTE: for simplicity, nzb_s_inner is used below, although strictly
1208!--          ----  speaking the following k-loop would have to be split up and
1209!--                rearranged according to the staggered grid.
1210             DO  k = nzb, nzt
1211                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j,i), 22 ) )
1212                ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +                  &
1213                                 u(k+1,j,i) - hom(k+1,1,1,sr) )
1214                vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +                  &
1215                                 v(k+1,j,i) - hom(k+1,1,2,sr) )
1216                pts = 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) +                 &
1217                                 pt(k+1,j,i) - hom(k+1,1,4,sr) )
1218
1219!--             Higher moments
1220                !$ACC ATOMIC
1221                sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 *        &
1222                                                    rmask(j,i,sr) * flag
1223                !$ACC ATOMIC
1224                sums_l(k,36,tn) = sums_l(k,36,tn) + pts**2 * w(k,j,i) *        &
1225                                                    rmask(j,i,sr) * flag
1226
1227!
1228!--             Salinity flux and density (density does not belong to here,
1229!--             but so far there is no other suitable place to calculate)
1230#ifndef _OPENACC
1231                IF ( ocean_mode )  THEN
1232                   IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
1233                      pts = 0.5_wp * ( sa(k,j,i)   - hom(k,1,23,sr) +          &
1234                                       sa(k+1,j,i) - hom(k+1,1,23,sr) )
1235                      sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) *     &
1236                                        rmask(j,i,sr) * flag
1237                   ENDIF
1238                   sums_l(k,64,tn) = sums_l(k,64,tn) + rho_ocean(k,j,i) *      &
1239                                                       rmask(j,i,sr) * flag
1240                   sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) *           &
1241                                                       rmask(j,i,sr) * flag
1242                ENDIF
1243
1244!
1245!--             Buoyancy flux, water flux, humidity flux, liquid water
1246!--             content, rain drop concentration and rain water content
1247                IF ( humidity )  THEN
1248                   IF ( bulk_cloud_model .OR. cloud_droplets )  THEN
1249                      pts = 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +         &
1250                                    vpt(k+1,j,i) - hom(k+1,1,44,sr) )
1251                      sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *     &
1252                                               rho_air_zw(k) *                 &
1253                                               heatflux_output_conversion(k) * &
1254                                                          rmask(j,i,sr) * flag
1255                      sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * rmask(j,i,sr) &
1256                                                                    * flag
1257
1258                      IF ( .NOT. cloud_droplets )  THEN
1259                         pts = 0.5_wp *                                        &
1260                              ( ( q(k,j,i) - ql(k,j,i) ) -                     &
1261                              hom(k,1,42,sr) +                                 &
1262                              ( q(k+1,j,i) - ql(k+1,j,i) ) -                   &
1263                              hom(k+1,1,42,sr) )
1264                         sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) *  &
1265                                             rho_air_zw(k) *                   &
1266                                             waterflux_output_conversion(k) *  &
1267                                                             rmask(j,i,sr)  *  &
1268                                                             flag
1269                         sums_l(k,75,tn) = sums_l(k,75,tn) + qc(k,j,i) *       &
1270                                                             rmask(j,i,sr) *   &
1271                                                             flag
1272                         sums_l(k,76,tn) = sums_l(k,76,tn) + prr(k,j,i) *      &
1273                                                             rmask(j,i,sr) *   &
1274                                                             flag
1275                         IF ( microphysics_morrison )  THEN
1276                            sums_l(k,123,tn) = sums_l(k,123,tn) + nc(k,j,i) *  &
1277                                                                rmask(j,i,sr) *&
1278                                                                flag
1279                         ENDIF
1280                         IF ( microphysics_seifert )  THEN
1281                            sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) *    &
1282                                                                rmask(j,i,sr) *&
1283                                                                flag 
1284                            sums_l(k,74,tn) = sums_l(k,74,tn) + qr(k,j,i) *    &
1285                                                                rmask(j,i,sr) *&
1286                                                                flag
1287                         ENDIF
1288                      ENDIF
1289
1290                   ELSE
1291                      IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
1292                         pts = 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +      &
1293                                          vpt(k+1,j,i) - hom(k+1,1,44,sr) )
1294                         sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *  &
1295                                              rho_air_zw(k) *                  &
1296                                              heatflux_output_conversion(k) *  &
1297                                                             rmask(j,i,sr)  *  &
1298                                                             flag
1299                      ELSE IF ( ws_scheme_sca .AND. sr == 0 )  THEN
1300                         sums_l(k,46,tn) = ( ( 1.0_wp + 0.61_wp *              & 
1301                                               hom(k,1,41,sr) ) *              &
1302                                             sums_l(k,17,tn) +                 &
1303                                             0.61_wp * hom(k,1,4,sr) *         &
1304                                             sums_l(k,49,tn)                   &
1305                                           ) * heatflux_output_conversion(k) * &
1306                                               flag
1307                      END IF
1308                   END IF
1309                ENDIF
1310!
1311!--             Passive scalar flux
1312                IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca                &
1313                     .OR. sr /= 0 ) )  THEN
1314                   pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,115,sr) +             &
1315                                    s(k+1,j,i) - hom(k+1,1,115,sr) )
1316                   sums_l(k,114,tn) = sums_l(k,114,tn) + pts * w(k,j,i) *      &
1317                                                       rmask(j,i,sr) * flag
1318                ENDIF
1319#endif
1320
1321!
1322!--             Energy flux w*e*
1323!--             has to be adjusted
1324                !$ACC ATOMIC
1325                sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5_wp *        &
1326                                             ( ust**2 + vst**2 + w(k,j,i)**2 ) &
1327                                           * rho_air_zw(k)                     &
1328                                           * momentumflux_output_conversion(k) &
1329                                           * rmask(j,i,sr) * flag
1330             ENDDO
1331          ENDDO
1332       ENDDO
1333       !$OMP END PARALLEL
1334
1335       !$ACC UPDATE &
1336       !$ACC HOST(sums_l(:,12,tn), sums_l(:,14,tn), sums_l(:,16,tn)) &
1337       !$ACC HOST(sums_l(:,35,tn), sums_l(:,36,tn), sums_l(:,37,tn))
1338!
1339!--    Treat land-surface quantities according to new wall model structure.
1340       IF ( land_surface )  THEN
1341          tn = 0
1342          !$OMP PARALLEL PRIVATE( i, j, m, tn )
1343          !$ tn = omp_get_thread_num()
1344          !$OMP DO
1345          DO  m = 1, surf_lsm_h%ns
1346             i = surf_lsm_h%i(m)
1347             j = surf_lsm_h%j(m)
1348       
1349             IF ( i >= nxl  .AND.  i <= nxr  .AND.                             &
1350                  j >= nys  .AND.  j <= nyn )  THEN
1351                sums_l(nzb,93,tn)  = sums_l(nzb,93,tn) + surf_lsm_h%ghf(m)
1352                sums_l(nzb,94,tn)  = sums_l(nzb,94,tn) + surf_lsm_h%qsws_liq(m)
1353                sums_l(nzb,95,tn)  = sums_l(nzb,95,tn) + surf_lsm_h%qsws_soil(m)
1354                sums_l(nzb,96,tn)  = sums_l(nzb,96,tn) + surf_lsm_h%qsws_veg(m)
1355                sums_l(nzb,97,tn)  = sums_l(nzb,97,tn) + surf_lsm_h%r_a(m)
1356                sums_l(nzb,98,tn) = sums_l(nzb,98,tn)+ surf_lsm_h%r_s(m)
1357             ENDIF
1358          ENDDO
1359          !$OMP END PARALLEL
1360
1361          tn = 0
1362          !$OMP PARALLEL PRIVATE( i, j, k, m, tn )
1363          !$ tn = omp_get_thread_num()
1364          !$OMP DO
1365          DO  m = 1, surf_lsm_h%ns
1366
1367             i = surf_lsm_h%i(m)           
1368             j = surf_lsm_h%j(m)
1369
1370             IF ( i >= nxl  .AND.  i <= nxr  .AND.                             &
1371                  j >= nys  .AND.  j <= nyn )  THEN
1372
1373                DO  k = nzb_soil, nzt_soil
1374                   sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil_h%var_2d(k,m)  &
1375                                      * rmask(j,i,sr)
1376                   sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil_h%var_2d(k,m)  &
1377                                      * rmask(j,i,sr)
1378                ENDDO
1379             ENDIF
1380          ENDDO
1381          !$OMP END PARALLEL
1382       ENDIF
1383!
1384!--    For speed optimization fluxes which have been computed in part directly
1385!--    inside the WS advection routines are treated seperatly
1386!--    Momentum fluxes first:
1387
1388       tn = 0
1389       !$OMP PARALLEL PRIVATE( i, j, k, tn, flag )
1390       !$ tn = omp_get_thread_num()
1391       IF ( .NOT. ws_scheme_mom .OR. sr /= 0  )  THEN
1392          !$OMP DO
1393          DO  i = nxl, nxr
1394             DO  j = nys, nyn
1395                DO  k = nzb, nzt
1396!
1397!--                Flag 23 is used to mask surface fluxes as well as model-top
1398!--                fluxes, which are added further below.
1399                   flag = MERGE( 1.0_wp, 0.0_wp,                               &
1400                                 BTEST( wall_flags_static_0(k,j,i), 23 ) ) *          &
1401                          MERGE( 1.0_wp, 0.0_wp,                               &
1402                                 BTEST( wall_flags_static_0(k,j,i), 9  ) )
1403
1404                   ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +               &
1405                                    u(k+1,j,i) - hom(k+1,1,1,sr) )
1406                   vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +               &
1407                                    v(k+1,j,i) - hom(k+1,1,2,sr) )
1408!
1409!--                Momentum flux w*u*
1410                   sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5_wp *                &
1411                                                     ( w(k,j,i-1) + w(k,j,i) ) &
1412                                           * rho_air_zw(k)                     &
1413                                           * momentumflux_output_conversion(k) &
1414                                                     * ust * rmask(j,i,sr)     &
1415                                                           * flag
1416!
1417!--                Momentum flux w*v*
1418                   sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5_wp *                &
1419                                                     ( w(k,j-1,i) + w(k,j,i) ) &
1420                                           * rho_air_zw(k)                     &
1421                                           * momentumflux_output_conversion(k) &
1422                                                     * vst * rmask(j,i,sr)     &
1423                                                           * flag
1424                ENDDO
1425             ENDDO
1426          ENDDO
1427
1428       ENDIF
1429       IF ( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
1430          !$OMP DO
1431          DO  i = nxl, nxr
1432             DO  j = nys, nyn
1433                DO  k = nzb, nzt
1434                   flag = MERGE( 1.0_wp, 0.0_wp,                               &
1435                                 BTEST( wall_flags_static_0(k,j,i), 23 ) ) *          &
1436                          MERGE( 1.0_wp, 0.0_wp,                               &
1437                                 BTEST( wall_flags_static_0(k,j,i), 9  ) )
1438!
1439!--                Vertical heat flux
1440                   sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5_wp *                &
1441                           ( pt(k,j,i)   - hom(k,1,4,sr) +                     &
1442                             pt(k+1,j,i) - hom(k+1,1,4,sr) )                   &
1443                           * rho_air_zw(k)                                     &
1444                           * heatflux_output_conversion(k)                     &
1445                           * w(k,j,i) * rmask(j,i,sr) * flag
1446                   IF ( humidity )  THEN
1447                      pts = 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +           &
1448                                      q(k+1,j,i) - hom(k+1,1,41,sr) )
1449                      sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *     &
1450                                       rho_air_zw(k) *                         &
1451                                       waterflux_output_conversion(k) *        &
1452                                       rmask(j,i,sr) * flag
1453                   ENDIF
1454                   IF ( passive_scalar )  THEN
1455                      pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,115,sr) +           &
1456                                      s(k+1,j,i) - hom(k+1,1,115,sr) )
1457                      sums_l(k,114,tn) = sums_l(k,114,tn) + pts * w(k,j,i) *    &
1458                                        rmask(j,i,sr) * flag
1459                   ENDIF
1460                ENDDO
1461             ENDDO
1462          ENDDO
1463
1464       ENDIF
1465
1466!
1467!--    Density at top follows Neumann condition
1468       IF ( ocean_mode )  THEN
1469          sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
1470          sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
1471       ENDIF
1472
1473!
1474!--    Divergence of vertical flux of resolved scale energy and pressure
1475!--    fluctuations as well as flux of pressure fluctuation itself (68).
1476!--    First calculate the products, then the divergence.
1477!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
1478       IF ( hom(nzb+1,2,55,0) /= 0.0_wp  .OR.  hom(nzb+1,2,68,0) /= 0.0_wp )   &
1479       THEN
1480          sums_ll = 0.0_wp  ! local array
1481
1482          !$OMP DO
1483          DO  i = nxl, nxr
1484             DO  j = nys, nyn
1485                DO  k = nzb+1, nzt
1486                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j,i), 0 ) )
1487
1488                   sums_ll(k,1) = sums_ll(k,1) + 0.5_wp * w(k,j,i) * (         &
1489                  ( 0.25_wp * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1) )  &
1490                            - 0.5_wp * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) ) )**2&
1491                + ( 0.25_wp * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i) )  &
1492                            - 0.5_wp * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) ) )**2&
1493                + w(k,j,i)**2                                        ) * flag
1494
1495                   sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i)             &
1496                                       * ( ( p(k,j,i) + p(k+1,j,i) )           &
1497                                         / momentumflux_output_conversion(k) ) &
1498                                       * flag
1499
1500                ENDDO
1501             ENDDO
1502          ENDDO
1503          sums_ll(0,1)     = 0.0_wp    ! because w is zero at the bottom
1504          sums_ll(nzt+1,1) = 0.0_wp
1505          sums_ll(0,2)     = 0.0_wp
1506          sums_ll(nzt+1,2) = 0.0_wp
1507
1508          DO  k = nzb+1, nzt
1509             sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
1510             sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
1511             sums_l(k,68,tn) = sums_ll(k,2)
1512          ENDDO
1513          sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
1514          sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
1515          sums_l(nzb,68,tn) = 0.0_wp    ! because w* = 0 at nzb
1516
1517       ENDIF
1518
1519!
1520!--    Divergence of vertical flux of SGS TKE and the flux itself (69)
1521       IF ( hom(nzb+1,2,57,0) /= 0.0_wp  .OR.  hom(nzb+1,2,69,0) /= 0.0_wp )   &
1522       THEN
1523          !$OMP DO
1524          DO  i = nxl, nxr
1525             DO  j = nys, nyn
1526                DO  k = nzb+1, nzt
1527
1528                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j,i), 0 ) )
1529
1530                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5_wp * (              &
1531                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
1532                 - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
1533                                                                ) * ddzw(k)    &
1534                                                                  * flag
1535
1536                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5_wp * (              &
1537                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
1538                                                                )  * flag
1539
1540                ENDDO
1541             ENDDO
1542          ENDDO
1543          sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
1544          sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
1545
1546       ENDIF
1547
1548!
1549!--    Horizontal heat fluxes (subgrid, resolved, total).
1550!--    Do it only, if profiles shall be plotted.
1551       IF ( hom(nzb+1,2,58,0) /= 0.0_wp ) THEN
1552
1553          !$OMP DO
1554          DO  i = nxl, nxr
1555             DO  j = nys, nyn
1556                DO  k = nzb+1, nzt
1557                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j,i), 0 ) )
1558!
1559!--                Subgrid horizontal heat fluxes u"pt", v"pt"
1560                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5_wp *                &
1561                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
1562                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
1563                                               * rho_air_zw(k)                 &
1564                                               * heatflux_output_conversion(k) &
1565                                                 * ddx * rmask(j,i,sr) * flag
1566                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
1567                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
1568                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
1569                                               * rho_air_zw(k)                 &
1570                                               * heatflux_output_conversion(k) &
1571                                                 * ddy * rmask(j,i,sr) * flag
1572!
1573!--                Resolved horizontal heat fluxes u*pt*, v*pt*
1574                   sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
1575                                                  ( u(k,j,i) - hom(k,1,1,sr) ) &
1576                                    * 0.5_wp * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
1577                                                 pt(k,j,i)   - hom(k,1,4,sr) ) &
1578                                               * heatflux_output_conversion(k) &
1579                                               * flag
1580                   pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
1581                                    pt(k,j,i)   - hom(k,1,4,sr) )
1582                   sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
1583                                                  ( v(k,j,i) - hom(k,1,2,sr) ) &
1584                                    * 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
1585                                                 pt(k,j,i)   - hom(k,1,4,sr) ) &
1586                                               * heatflux_output_conversion(k) &
1587                                               * flag
1588                ENDDO
1589             ENDDO
1590          ENDDO
1591!
1592!--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
1593          sums_l(nzb,58,tn) = 0.0_wp
1594          sums_l(nzb,59,tn) = 0.0_wp
1595          sums_l(nzb,60,tn) = 0.0_wp
1596          sums_l(nzb,61,tn) = 0.0_wp
1597          sums_l(nzb,62,tn) = 0.0_wp
1598          sums_l(nzb,63,tn) = 0.0_wp
1599
1600       ENDIF
1601       !$OMP END PARALLEL
1602
1603!
1604!--    Collect current large scale advection and subsidence tendencies for
1605!--    data output
1606       IF ( large_scale_forcing  .AND.  ( simulated_time > 0.0_wp ) )  THEN
1607!
1608!--       Interpolation in time of LSF_DATA
1609          nt = 1
1610          DO WHILE ( simulated_time - dt_3d > time_vert(nt) )
1611             nt = nt + 1
1612          ENDDO
1613          IF ( simulated_time - dt_3d /= time_vert(nt) )  THEN
1614            nt = nt - 1
1615          ENDIF
1616
1617          fac = ( simulated_time - dt_3d - time_vert(nt) )                     &
1618                / ( time_vert(nt+1)-time_vert(nt) )
1619
1620
1621          DO  k = nzb, nzt
1622             sums_ls_l(k,0) = td_lsa_lpt(k,nt)                                 &
1623                              + fac * ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )
1624             sums_ls_l(k,1) = td_lsa_q(k,nt)                                   &
1625                              + fac * ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) )
1626          ENDDO
1627
1628          sums_ls_l(nzt+1,0) = sums_ls_l(nzt,0)
1629          sums_ls_l(nzt+1,1) = sums_ls_l(nzt,1)
1630
1631          IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
1632
1633             DO  k = nzb, nzt
1634                sums_ls_l(k,2) = td_sub_lpt(k,nt) + fac *                      &
1635                                 ( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )
1636                sums_ls_l(k,3) = td_sub_q(k,nt) + fac *                        &
1637                                 ( td_sub_q(k,nt+1) - td_sub_q(k,nt) )
1638             ENDDO
1639
1640             sums_ls_l(nzt+1,2) = sums_ls_l(nzt,2)
1641             sums_ls_l(nzt+1,3) = sums_ls_l(nzt,3)
1642
1643          ENDIF
1644
1645       ENDIF
1646
1647       tn = 0
1648       !$OMP PARALLEL PRIVATE( i, j, k, tn )
1649       !$ tn = omp_get_thread_num()       
1650       IF ( radiation .AND. radiation_scheme == 'rrtmg' )  THEN
1651          !$OMP DO
1652          DO  i = nxl, nxr
1653             DO  j =  nys, nyn
1654                DO  k = nzb+1, nzt+1
1655                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_static_0(k,j,i), 0 ) )
1656
1657                   sums_l(k,100,tn)  = sums_l(k,100,tn)  + rad_lw_in(k,j,i)    &
1658                                       * rmask(j,i,sr) * flag
1659                   sums_l(k,101,tn)  = sums_l(k,101,tn)  + rad_lw_out(k,j,i)   &
1660                                       * rmask(j,i,sr) * flag
1661                   sums_l(k,102,tn)  = sums_l(k,102,tn)  + rad_sw_in(k,j,i)    &
1662                                       * rmask(j,i,sr) * flag
1663                   sums_l(k,103,tn)  = sums_l(k,103,tn)  + rad_sw_out(k,j,i)   &
1664                                       * rmask(j,i,sr) * flag
1665                   sums_l(k,104,tn)  = sums_l(k,104,tn)  + rad_lw_cs_hr(k,j,i) &
1666                                       * rmask(j,i,sr) * flag
1667                   sums_l(k,105,tn)  = sums_l(k,105,tn)  + rad_lw_hr(k,j,i)    &
1668                                       * rmask(j,i,sr) * flag
1669                   sums_l(k,106,tn)  = sums_l(k,106,tn)  + rad_sw_cs_hr(k,j,i) &
1670                                       * rmask(j,i,sr) * flag
1671                   sums_l(k,107,tn)  = sums_l(k,107,tn)  + rad_sw_hr(k,j,i)    &
1672                                       * rmask(j,i,sr) * flag
1673                ENDDO
1674             ENDDO
1675          ENDDO
1676       ENDIF
1677
1678!
1679!--    Calculate the profiles for all other modules
1680       CALL module_interface_statistics( 'profiles', sr, tn, dots_max )
1681       !$OMP END PARALLEL
1682
1683!
1684!--    Summation of thread sums
1685       IF ( threads_per_task > 1 )  THEN
1686          DO  i = 1, threads_per_task-1
1687             sums_l(:,3,0)          = sums_l(:,3,0) + sums_l(:,3,i)
1688             sums_l(:,4:40,0)       = sums_l(:,4:40,0) + sums_l(:,4:40,i)
1689             sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
1690                                      sums_l(:,45:pr_palm,i)
1691             IF ( max_pr_user > 0 )  THEN
1692                sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
1693                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
1694                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
1695             ENDIF
1696
1697             IF ( air_chemistry )  THEN
1698                IF ( max_pr_cs > 0 )  THEN                                 
1699                     sums_l(:,pr_palm+max_pr_user+1:pr_palm + max_pr_user+ max_pr_cs,0) =          &
1700                               sums_l(:,pr_palm+max_pr_user+1:pr_palm + max_pr_user+max_pr_cs,0) + &
1701                               sums_l(:,pr_palm+max_pr_user+1:pr_palm + max_pr_user+max_pr_cs,i)
1702
1703                ENDIF
1704             ENDIF
1705             IF ( salsa )  THEN
1706                IF ( max_pr_cs > 0 )  THEN
1707                   sums_l(:,pr_palm+max_pr_user+max_pr_cs+1:pr_palm+max_pr_user+max_pr_cs+max_pr_salsa,0) =    &
1708                      sums_l(:,pr_palm+max_pr_user+max_pr_cs+1:pr_palm+max_pr_user+max_pr_cs+max_pr_salsa,0) + &
1709                      sums_l(:,pr_palm+max_pr_user+max_pr_cs+1:pr_palm+max_pr_user+max_pr_cs+max_pr_salsa,i)
1710
1711                ENDIF
1712             ENDIF
1713          ENDDO
1714       ENDIF
1715
1716#if defined( __parallel )
1717
1718!
1719!--    Compute total sum from local sums
1720       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1721       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL,   &
1722                           MPI_SUM, comm2d, ierr )
1723       IF ( large_scale_forcing )  THEN
1724          CALL MPI_ALLREDUCE( sums_ls_l(nzb,2), sums(nzb,83), ngp_sums_ls,     &
1725                              MPI_REAL, MPI_SUM, comm2d, ierr )
1726       ENDIF
1727
1728       IF ( air_chemistry  .AND.  max_pr_cs > 0 )  THEN
1729          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1730             DO  i = 1, max_pr_cs
1731                CALL MPI_ALLREDUCE( sums_l(nzb,pr_palm+max_pr_user+i,0),       &
1732                                    sums(nzb,pr_palm+max_pr_user+i),           &
1733                                    nzt+2-nzb, MPI_REAL, MPI_SUM, comm2d, ierr )
1734             ENDDO
1735       ENDIF
1736
1737       IF ( salsa  .AND.  max_pr_salsa > 0 )  THEN
1738          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1739             DO  i = 1, max_pr_salsa
1740                CALL MPI_ALLREDUCE( sums_l(nzb,pr_palm+max_pr_user+max_pr_cs+i,0),                 &
1741                                    sums(nzb,pr_palm+max_pr_user+max_pr_user+i),                   &
1742                                    nzt+2-nzb, MPI_REAL, MPI_SUM, comm2d, ierr )
1743             ENDDO
1744       ENDIF
1745
1746#else
1747       sums = sums_l(:,:,0)
1748       IF ( large_scale_forcing )  THEN
1749          sums(:,81:88) = sums_ls_l
1750       ENDIF
1751#endif
1752
1753!
1754!--    Final values are obtained by division by the total number of grid points
1755!--    used for summation. After that store profiles.
1756!--    Check, if statistical regions do contain at least one grid point at the
1757!--    respective k-level, otherwise division by zero will lead to undefined
1758!--    values, which may cause e.g. problems with NetCDF output
1759!--    Profiles:
1760       DO  k = nzb, nzt+1
1761          sums(k,3)             = sums(k,3)             / ngp_2dh(sr)
1762          sums(k,12:22)         = sums(k,12:22)         / ngp_2dh(sr)
1763          sums(k,30:32)         = sums(k,30:32)         / ngp_2dh(sr)
1764          sums(k,35:39)         = sums(k,35:39)         / ngp_2dh(sr)
1765          sums(k,45:53)         = sums(k,45:53)         / ngp_2dh(sr)
1766          sums(k,55:63)         = sums(k,55:63)         / ngp_2dh(sr)
1767          sums(k,81:88)         = sums(k,81:88)         / ngp_2dh(sr)
1768          sums(k,89:112)        = sums(k,89:112)        / ngp_2dh(sr)
1769          sums(k,114)           = sums(k,114)           / ngp_2dh(sr)
1770          sums(k,117)           = sums(k,117)           / ngp_2dh(sr)
1771          IF ( ngp_2dh_s_inner(k,sr) /= 0 )  THEN
1772             sums(k,8:11)          = sums(k,8:11)          / ngp_2dh_s_inner(k,sr)
1773             sums(k,23:29)         = sums(k,23:29)         / ngp_2dh_s_inner(k,sr)
1774             sums(k,33:34)         = sums(k,33:34)         / ngp_2dh_s_inner(k,sr)
1775             sums(k,40)            = sums(k,40)            / ngp_2dh_s_inner(k,sr)
1776             sums(k,54)            = sums(k,54)            / ngp_2dh_s_inner(k,sr)
1777             sums(k,64)            = sums(k,64)            / ngp_2dh_s_inner(k,sr)
1778             sums(k,70:80)         = sums(k,70:80)         / ngp_2dh_s_inner(k,sr)
1779             sums(k,116)           = sums(k,116)           / ngp_2dh_s_inner(k,sr)
1780             sums(k,118:pr_palm-2) = sums(k,118:pr_palm-2) / ngp_2dh_s_inner(k,sr)
1781             sums(k,123)           = sums(k,123) * ngp_2dh_s_inner(k,sr)  / ngp_2dh(sr)
1782          ENDIF
1783       ENDDO
1784
1785!--    u* and so on
1786!--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
1787!--    size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
1788!--    above the topography, they are being divided by ngp_2dh(sr)
1789       sums(nzb:nzb+3,pr_palm)    = sums(nzb:nzb+3,pr_palm)    / &
1790                                    ngp_2dh(sr)
1791       sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)       / &    ! qs
1792                                    ngp_2dh(sr)
1793       sums(nzb+13,pr_palm)       = sums(nzb+13,pr_palm)       / &    ! ss
1794                                    ngp_2dh(sr)
1795       sums(nzb+14,pr_palm)       = sums(nzb+14,pr_palm)       / &    ! surface temperature
1796                                    ngp_2dh(sr)
1797!--    eges, e*
1798       sums(nzb+4:nzb+5,pr_palm)  = sums(nzb+4:nzb+5,pr_palm)  / &
1799                                    ngp_3d(sr)
1800!--    Old and new divergence
1801       sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
1802                                    ngp_3d_inner(sr)
1803
1804!--    User-defined profiles
1805       IF ( max_pr_user > 0 )  THEN
1806          DO  k = nzb, nzt+1
1807             sums(k,pr_palm+1:pr_palm+max_pr_user) = &
1808                                    sums(k,pr_palm+1:pr_palm+max_pr_user) / &
1809                                    ngp_2dh_s_inner(k,sr)
1810          ENDDO
1811       ENDIF
1812
1813       IF ( air_chemistry ) THEN
1814          IF ( max_pr_cs > 0 )  THEN                 
1815             DO k = nzb, nzt+1
1816                sums(k, pr_palm+1:pr_palm+max_pr_user+max_pr_cs) = &
1817                                 sums(k, pr_palm+1:pr_palm+max_pr_user+max_pr_cs) / &
1818                                 ngp_2dh_s_inner(k,sr)
1819             ENDDO
1820          ENDIF 
1821       ENDIF
1822
1823       IF ( salsa ) THEN
1824          IF ( max_pr_salsa > 0 )  THEN
1825             DO k = nzb, nzt+1
1826                sums(k,pr_palm+max_pr_user+max_pr_cs+1:pr_palm+max_pr_user+max_pr_cs+max_pr_salsa) = &
1827                  sums(k,pr_palm+max_pr_user+max_pr_cs+1:pr_palm+max_pr_user+max_pr_cs+max_pr_salsa) &
1828                  / ngp_2dh_s_inner(k,sr)
1829             ENDDO
1830          ENDIF 
1831       ENDIF
1832
1833!
1834!--    Collect horizontal average in hom.
1835!--    Compute deduced averages (e.g. total heat flux)
1836       hom(:,1,3,sr)  = sums(:,3)      ! w
1837       hom(:,1,8,sr)  = sums(:,8)      ! e     profiles 5-7 are initial profiles
1838       hom(:,1,9,sr)  = sums(:,9)      ! km
1839       hom(:,1,10,sr) = sums(:,10)     ! kh
1840       hom(:,1,11,sr) = sums(:,11)     ! l
1841       hom(:,1,12,sr) = sums(:,12)     ! w"u"
1842       hom(:,1,13,sr) = sums(:,13)     ! w*u*
1843       hom(:,1,14,sr) = sums(:,14)     ! w"v"
1844       hom(:,1,15,sr) = sums(:,15)     ! w*v*
1845       hom(:,1,16,sr) = sums(:,16)     ! w"pt"
1846       hom(:,1,17,sr) = sums(:,17)     ! w*pt*
1847       hom(:,1,18,sr) = sums(:,16) + sums(:,17)    ! wpt
1848       hom(:,1,19,sr) = sums(:,12) + sums(:,13)    ! wu
1849       hom(:,1,20,sr) = sums(:,14) + sums(:,15)    ! wv
1850       hom(:,1,21,sr) = sums(:,21)     ! w*pt*BC
1851       hom(:,1,22,sr) = sums(:,16) + sums(:,21)    ! wptBC
1852                                       ! profile 24 is initial profile (sa)
1853                                       ! profiles 25-29 left empty for initial
1854                                       ! profiles
1855       hom(:,1,30,sr) = sums(:,30)     ! u*2
1856       hom(:,1,31,sr) = sums(:,31)     ! v*2
1857       hom(:,1,32,sr) = sums(:,32)     ! w*2
1858       hom(:,1,33,sr) = sums(:,33)     ! pt*2
1859       hom(:,1,34,sr) = sums(:,34)     ! e*
1860       hom(:,1,35,sr) = sums(:,35)     ! w*2pt*
1861       hom(:,1,36,sr) = sums(:,36)     ! w*pt*2
1862       hom(:,1,37,sr) = sums(:,37)     ! w*e*
1863       hom(:,1,38,sr) = sums(:,38)     ! w*3
1864       hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20_wp )**1.5_wp   ! Sw
1865       hom(:,1,40,sr) = sums(:,40)     ! p
1866       hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
1867       hom(:,1,46,sr) = sums(:,46)     ! w*vpt*       
1868       hom(:,1,47,sr) = sums(:,45) + sums(:,46)    ! wvpt
1869       hom(:,1,48,sr) = sums(:,48)     ! w"q" (w"qv")
1870       hom(:,1,49,sr) = sums(:,49)     ! w*q* (w*qv*)
1871       hom(:,1,50,sr) = sums(:,48) + sums(:,49)    ! wq (wqv)
1872       hom(:,1,51,sr) = sums(:,51)     ! w"qv"
1873       hom(:,1,52,sr) = sums(:,52)     ! w*qv*       
1874       hom(:,1,53,sr) = sums(:,52) + sums(:,51)    ! wq (wqv)
1875       hom(:,1,54,sr) = sums(:,54)     ! ql
1876       hom(:,1,55,sr) = sums(:,55)     ! w*u*u*/dz
1877       hom(:,1,56,sr) = sums(:,56)     ! w*p*/dz
1878       hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho_ocean )/dz
1879       hom(:,1,58,sr) = sums(:,58)     ! u"pt"
1880       hom(:,1,59,sr) = sums(:,59)     ! u*pt*
1881       hom(:,1,60,sr) = sums(:,58) + sums(:,59)    ! upt_t
1882       hom(:,1,61,sr) = sums(:,61)     ! v"pt"
1883       hom(:,1,62,sr) = sums(:,62)     ! v*pt*
1884       hom(:,1,63,sr) = sums(:,61) + sums(:,62)    ! vpt_t
1885       hom(:,1,64,sr) = sums(:,64)     ! rho_ocean
1886       hom(:,1,65,sr) = sums(:,65)     ! w"sa"
1887       hom(:,1,66,sr) = sums(:,66)     ! w*sa*
1888       hom(:,1,67,sr) = sums(:,65) + sums(:,66)    ! wsa
1889       hom(:,1,68,sr) = sums(:,68)     ! w*p*
1890       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho_ocean
1891       hom(:,1,70,sr) = sums(:,70)     ! q*2
1892       hom(:,1,71,sr) = sums(:,71)     ! prho
1893       hom(:,1,72,sr) = hyp * 1E-2_wp  ! hyp in hPa
1894       hom(:,1,123,sr) = sums(:,123)   ! nc
1895       hom(:,1,73,sr) = sums(:,73)     ! nr
1896       hom(:,1,74,sr) = sums(:,74)     ! qr
1897       hom(:,1,75,sr) = sums(:,75)     ! qc
1898       hom(:,1,76,sr) = sums(:,76)     ! prr (precipitation rate)
1899                                       ! 77 is initial density profile
1900       hom(:,1,78,sr) = ug             ! ug
1901       hom(:,1,79,sr) = vg             ! vg
1902       hom(:,1,80,sr) = w_subs         ! w_subs
1903
1904       IF ( large_scale_forcing )  THEN
1905          hom(:,1,81,sr) = sums_ls_l(:,0)          ! td_lsa_lpt
1906          hom(:,1,82,sr) = sums_ls_l(:,1)          ! td_lsa_q
1907          IF ( use_subsidence_tendencies )  THEN
1908             hom(:,1,83,sr) = sums_ls_l(:,2)       ! td_sub_lpt
1909             hom(:,1,84,sr) = sums_ls_l(:,3)       ! td_sub_q
1910          ELSE
1911             hom(:,1,83,sr) = sums(:,83)           ! td_sub_lpt
1912             hom(:,1,84,sr) = sums(:,84)           ! td_sub_q
1913          ENDIF
1914          hom(:,1,85,sr) = sums(:,85)              ! td_nud_lpt
1915          hom(:,1,86,sr) = sums(:,86)              ! td_nud_q
1916          hom(:,1,87,sr) = sums(:,87)              ! td_nud_u
1917          hom(:,1,88,sr) = sums(:,88)              ! td_nud_v
1918       ENDIF
1919
1920       IF ( land_surface )  THEN
1921          hom(:,1,89,sr) = sums(:,89)              ! t_soil
1922                                                   ! 90 is initial t_soil profile
1923          hom(:,1,91,sr) = sums(:,91)              ! m_soil
1924                                                   ! 92 is initial m_soil profile
1925          hom(:,1,93,sr)  = sums(:,93)             ! ghf
1926          hom(:,1,94,sr)  = sums(:,94)             ! qsws_liq
1927          hom(:,1,95,sr)  = sums(:,95)             ! qsws_soil
1928          hom(:,1,96,sr)  = sums(:,96)             ! qsws_veg
1929          hom(:,1,97,sr)  = sums(:,97)             ! r_a
1930          hom(:,1,98,sr) = sums(:,98)              ! r_s
1931
1932       ENDIF
1933
1934       IF ( radiation )  THEN
1935          hom(:,1,99,sr) = sums(:,99)            ! rad_net
1936          hom(:,1,100,sr) = sums(:,100)            ! rad_lw_in
1937          hom(:,1,101,sr) = sums(:,101)            ! rad_lw_out
1938          hom(:,1,102,sr) = sums(:,102)            ! rad_sw_in
1939          hom(:,1,103,sr) = sums(:,103)            ! rad_sw_out
1940
1941          IF ( radiation_scheme == 'rrtmg' )  THEN
1942             hom(:,1,104,sr) = sums(:,104)            ! rad_lw_cs_hr
1943             hom(:,1,105,sr) = sums(:,105)            ! rad_lw_hr
1944             hom(:,1,106,sr) = sums(:,106)            ! rad_sw_cs_hr
1945             hom(:,1,107,sr) = sums(:,107)            ! rad_sw_hr
1946
1947             hom(:,1,108,sr) = sums(:,108)            ! rrtm_aldif
1948             hom(:,1,109,sr) = sums(:,109)            ! rrtm_aldir
1949             hom(:,1,110,sr) = sums(:,110)            ! rrtm_asdif
1950             hom(:,1,111,sr) = sums(:,111)            ! rrtm_asdir
1951          ENDIF
1952       ENDIF
1953
1954       hom(:,1,112,sr) = sums(:,112)            !: L
1955
1956       IF ( passive_scalar )  THEN
1957          hom(:,1,117,sr) = sums(:,117)     ! w"s"
1958          hom(:,1,114,sr) = sums(:,114)     ! w*s*
1959          hom(:,1,118,sr) = sums(:,117) + sums(:,114)    ! ws
1960          hom(:,1,116,sr) = sums(:,116)     ! s*2
1961       ENDIF
1962
1963       hom(:,1,119,sr) = rho_air       ! rho_air in Kg/m^3
1964       hom(:,1,120,sr) = rho_air_zw    ! rho_air_zw in Kg/m^3
1965
1966       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
1967                                       ! u*, w'u', w'v', t* (in last profile)
1968
1969       IF ( max_pr_user > 0 )  THEN    ! user-defined profiles
1970          hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
1971                               sums(:,pr_palm+1:pr_palm+max_pr_user)
1972       ENDIF
1973
1974       IF ( air_chemistry )  THEN
1975          IF ( max_pr_cs > 0 )  THEN    ! chem_spcs profiles     
1976             hom(:, 1, pr_palm+max_pr_user+1:pr_palm + max_pr_user+max_pr_cs, sr) = &
1977                               sums(:, pr_palm+max_pr_user+1:pr_palm+max_pr_user+max_pr_cs)
1978          ENDIF
1979       ENDIF
1980
1981       IF ( salsa )  THEN
1982          IF ( max_pr_salsa > 0 )  THEN    ! salsa profiles
1983             hom(:,1,pr_palm+max_pr_user+max_pr_cs+1:pr_palm+max_pr_user+max_pr_cs+max_pr_salsa, sr) = &
1984                  sums(:,pr_palm+max_pr_user+max_pr_cs+1:pr_palm+max_pr_user+max_pr_cs+max_pr_salsa)
1985          ENDIF
1986       ENDIF
1987!
1988!--    Determine the boundary layer height using two different schemes.
1989!--    First scheme: Starting from the Earth's (Ocean's) surface, look for the
1990!--    first relative minimum (maximum) of the total heat flux.
1991!--    The corresponding height is assumed as the boundary layer height, if it
1992!--    is less than 1.5 times the height where the heat flux becomes negative
1993!--    (positive) for the first time. Attention: the resolved vertical sensible
1994!--    heat flux (hom(:,1,17,sr) = w*pt*) is not known at the beginning because
1995!--    the calculation happens in advec_s_ws which is called after
1996!--    flow_statistics. Therefore z_i is directly taken from restart data at
1997!--    the beginning of restart runs. 
1998       IF ( TRIM( initializing_actions ) /= 'read_restart_data' .OR.           &
1999            simulated_time_at_begin /= simulated_time ) THEN
2000
2001          z_i(1) = 0.0_wp
2002          first = .TRUE.
2003
2004          IF ( ocean_mode )  THEN
2005             DO  k = nzt, nzb+1, -1
2006                IF ( first  .AND.  hom(k,1,18,sr) < -1.0E-8_wp )  THEN
2007                   first = .FALSE.
2008                   height = zw(k)
2009                ENDIF
2010                IF ( hom(k,1,18,sr) < -1.0E-8_wp  .AND.                        &
2011                     hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
2012                   IF ( zw(k) < 1.5_wp * height )  THEN
2013                      z_i(1) = zw(k)
2014                   ELSE
2015                      z_i(1) = height
2016                   ENDIF
2017                   EXIT
2018                ENDIF
2019             ENDDO
2020          ELSE
2021             DO  k = nzb, nzt-1
2022                IF ( first  .AND.  hom(k,1,18,sr) < -1.0E-8_wp )  THEN
2023                   first = .FALSE.
2024                   height = zw(k)
2025                ENDIF
2026                IF ( hom(k,1,18,sr) < -1.0E-8_wp  .AND.                        &
2027                     hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
2028                   IF ( zw(k) < 1.5_wp * height )  THEN
2029                      z_i(1) = zw(k)
2030                   ELSE
2031                      z_i(1) = height
2032                   ENDIF
2033                   EXIT
2034                ENDIF
2035             ENDDO
2036          ENDIF
2037
2038!
2039!--       Second scheme: Gradient scheme from Sullivan et al. (1998), modified
2040!--       by Uhlenbrock(2006). The boundary layer height is the height with the
2041!--       maximal local temperature gradient: starting from the second (the
2042!--       last but one) vertical gridpoint, the local gradient must be at least
2043!--       0.2K/100m and greater than the next four gradients.
2044!--       WARNING: The threshold value of 0.2K/100m must be adjusted for the
2045!--       ocean case!
2046          z_i(2) = 0.0_wp
2047          DO  k = nzb+1, nzt+1
2048             dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
2049          ENDDO
2050          dptdz_threshold = 0.2_wp / 100.0_wp
2051
2052          IF ( ocean_mode )  THEN
2053             DO  k = nzt+1, nzb+5, -1
2054                IF ( dptdz(k) > dptdz_threshold  .AND.                         &
2055                     dptdz(k) > dptdz(k-1)  .AND.  dptdz(k) > dptdz(k-2)  .AND.&
2056                     dptdz(k) > dptdz(k-3)  .AND.  dptdz(k) > dptdz(k-4) )  THEN
2057                   z_i(2) = zw(k-1)
2058                   EXIT
2059                ENDIF
2060             ENDDO
2061          ELSE
2062             DO  k = nzb+1, nzt-3
2063                IF ( dptdz(k) > dptdz_threshold  .AND.                         &
2064                     dptdz(k) > dptdz(k+1)  .AND.  dptdz(k) > dptdz(k+2)  .AND.&
2065                     dptdz(k) > dptdz(k+3)  .AND.  dptdz(k) > dptdz(k+4) )  THEN
2066                   z_i(2) = zw(k-1)
2067                   EXIT
2068                ENDIF
2069             ENDDO
2070          ENDIF
2071
2072       ENDIF
2073
2074       hom(nzb+6,1,pr_palm,sr) = z_i(1)
2075       hom(nzb+7,1,pr_palm,sr) = z_i(2)
2076
2077!
2078!--    Determine vertical index which is nearest to the mean surface level
2079!--    height of the respective statistic region
2080       DO  k = nzb, nzt
2081          IF ( zw(k) >= mean_surface_level_height(sr) )  THEN
2082             k_surface_level = k
2083             EXIT
2084          ENDIF
2085       ENDDO
2086
2087!
2088!--    Computation of both the characteristic vertical velocity and
2089!--    the characteristic convective boundary layer temperature.
2090!--    The inversion height entering into the equation is defined with respect
2091!--    to the mean surface level height of the respective statistic region.
2092!--    The horizontal average at surface level index + 1 is input for the
2093!--    average temperature.
2094       IF ( hom(k_surface_level,1,18,sr) > 1.0E-8_wp  .AND.  z_i(1) /= 0.0_wp )&
2095       THEN
2096          hom(nzb+8,1,pr_palm,sr) =                                            &
2097             ( g / hom(k_surface_level+1,1,4,sr) *                             &
2098             ( hom(k_surface_level,1,18,sr) /                                  &
2099             ( heatflux_output_conversion(nzb) * rho_air(nzb) ) )              &
2100             * ABS( z_i(1) - mean_surface_level_height(sr) ) )**0.333333333_wp
2101       ELSE
2102          hom(nzb+8,1,pr_palm,sr)  = 0.0_wp
2103       ENDIF
2104
2105!
2106!--    Collect the time series quantities. Please note, timeseries quantities
2107!--    which are collected from horizontally averaged profiles, e.g. wpt
2108!--    or pt(zp), are treated specially. In case of elevated model surfaces,
2109!--    index nzb+1 might be within topography and data will be zero. Therefore,
2110!--    take value for the first atmosphere index, which is topo_min_level+1.
2111       ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)        ! E
2112       ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)        ! E*
2113       ts_value(3,sr) = dt_3d
2114       ts_value(4,sr) = hom(nzb,1,pr_palm,sr)          ! u*
2115       ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr)        ! th*
2116       ts_value(6,sr) = u_max
2117       ts_value(7,sr) = v_max
2118       ts_value(8,sr) = w_max
2119       ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr)       ! new divergence
2120       ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr)       ! old Divergence
2121       ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr)       ! z_i(1)
2122       ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr)       ! z_i(2)
2123       ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr)       ! w*
2124       ts_value(14,sr) = hom(nzb,1,16,sr)              ! w'pt'   at k=0
2125       ts_value(15,sr) = hom(topo_min_level+1,1,16,sr) ! w'pt'   at k=1
2126       ts_value(16,sr) = hom(topo_min_level+1,1,18,sr) ! wpt     at k=1
2127       ts_value(17,sr) = hom(nzb+14,1,pr_palm,sr)      ! pt(0)
2128       ts_value(18,sr) = hom(topo_min_level+1,1,4,sr)  ! pt(zp)
2129       ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr)       ! u'w'    at k=0
2130       ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)       ! v'w'    at k=0
2131       ts_value(21,sr) = hom(nzb,1,48,sr)              ! w"q"    at k=0
2132
2133       IF ( .NOT. neutral )  THEN
2134          ts_value(22,sr) = hom(nzb,1,112,sr)          ! L
2135       ELSE
2136          ts_value(22,sr) = 1.0E10_wp
2137       ENDIF
2138
2139       ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
2140
2141       IF ( passive_scalar )  THEN
2142          ts_value(24,sr) = hom(nzb+13,1,117,sr)       ! w"s" ( to do ! )
2143          ts_value(25,sr) = hom(nzb+13,1,pr_palm,sr)   ! s*
2144       ENDIF
2145
2146!
2147!--    Collect land surface model timeseries
2148       IF ( land_surface )  THEN
2149          ts_value(dots_soil  ,sr) = hom(nzb,1,93,sr)           ! ghf
2150          ts_value(dots_soil+1,sr) = hom(nzb,1,94,sr)           ! qsws_liq
2151          ts_value(dots_soil+2,sr) = hom(nzb,1,95,sr)           ! qsws_soil
2152          ts_value(dots_soil+3,sr) = hom(nzb,1,96,sr)           ! qsws_veg
2153          ts_value(dots_soil+4,sr) = hom(nzb,1,97,sr)           ! r_a
2154          ts_value(dots_soil+5,sr) = hom(nzb,1,98,sr)           ! r_s
2155       ENDIF
2156!
2157!--    Collect radiation model timeseries
2158       IF ( radiation )  THEN
2159          ts_value(dots_rad,sr)   = hom(nzb,1,99,sr)           ! rad_net
2160          ts_value(dots_rad+1,sr) = hom(nzb,1,100,sr)          ! rad_lw_in
2161          ts_value(dots_rad+2,sr) = hom(nzb,1,101,sr)          ! rad_lw_out
2162          ts_value(dots_rad+3,sr) = hom(nzb,1,102,sr)          ! rad_sw_in
2163          ts_value(dots_rad+4,sr) = hom(nzb,1,103,sr)          ! rad_sw_out
2164
2165          IF ( radiation_scheme == 'rrtmg' )  THEN
2166             ts_value(dots_rad+5,sr) = hom(nzb,1,108,sr)          ! rrtm_aldif
2167             ts_value(dots_rad+6,sr) = hom(nzb,1,109,sr)          ! rrtm_aldir
2168             ts_value(dots_rad+7,sr) = hom(nzb,1,110,sr)          ! rrtm_asdif
2169             ts_value(dots_rad+8,sr) = hom(nzb,1,111,sr)          ! rrtm_asdir
2170          ENDIF
2171
2172       ENDIF
2173
2174!
2175!--    Calculate additional statistics provided by other modules
2176       CALL module_interface_statistics( 'time_series', sr, 0, dots_max )
2177
2178    ENDDO    ! loop of the subregions
2179
2180!
2181!-- If required, sum up horizontal averages for subsequent time averaging.
2182!-- Do not sum, if flow statistics is called before the first initial time step.
2183    IF ( do_sum  .AND.  simulated_time /= 0.0_wp )  THEN
2184       IF ( average_count_pr == 0 )  hom_sum = 0.0_wp
2185       hom_sum = hom_sum + hom(:,1,:,:)
2186       average_count_pr = average_count_pr + 1
2187       do_sum = .FALSE.
2188    ENDIF
2189
2190!
2191!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
2192!-- This flag is reset after each time step in time_integration.
2193    flow_statistics_called = .TRUE.
2194
2195    CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
2196
2197
2198 END SUBROUTINE flow_statistics
Note: See TracBrowser for help on using the repository browser.