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

Last change on this file since 3651 was 3651, checked in by suehring, 2 years ago

Bugfix, restore OMP END PARALLEL block (accidantly remove in -r 3637)

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