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

Last change on this file since 3254 was 3241, checked in by raasch, 6 years ago

various changes to avoid compiler warnings (mainly removal of unused variables)

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