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

Last change on this file since 3534 was 3458, checked in by kanani, 5 years ago

Reintegrated fixes/changes from branch chemistry

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