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

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

Modularization of all bulk cloud physics code components

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