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

Last change on this file since 3452 was 3452, checked in by schwenkel, 5 years ago

Bugfix for profiles output

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