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

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

Bugfix in output conversion of resolved-scale momentum fluxes in case of PW advection scheme

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