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

Last change on this file since 2233 was 2233, checked in by suehring, 4 years ago

last commit documented

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