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

Last change on this file since 2894 was 2817, checked in by knoop, 6 years ago

Preliminary gust module interface implemented

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