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

Last change on this file since 2753 was 2753, checked in by suehring, 3 years ago

Tile approach for spectral albedo implemented.

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