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

Last change on this file since 2270 was 2270, checked in by maronga, 4 years ago

major revisions in land surface model

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