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

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