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

Last change on this file since 3371 was 3298, checked in by kanani, 6 years ago

Merge chemistry branch at r3297 to trunk

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