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

Last change on this file since 2232 was 2232, checked in by suehring, 4 years ago

Adjustments according new topography and surface-modelling concept implemented

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