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

Last change on this file since 3650 was 3637, checked in by knoop, 5 years ago

M Makefile

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