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

Last change on this file since 2968 was 2968, checked in by suehring, 6 years ago

Bugfixes in initalization of land-surface model as well as timeseries data output in case of elevated model surfaces.

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