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

Last change on this file since 3294 was 3294, checked in by raasch, 3 years ago

modularization of the ocean code

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