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

Last change on this file since 3022 was 3004, checked in by Giersch, 6 years ago

precipitation_rate removed, further allocation checks for data output of averaged quantities implemented, double CALL of flow_statistics at the beginning of time_integration removed, further minor bugfixes, comments added

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