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

Last change on this file since 3452 was 3452, checked in by schwenkel, 3 years ago

Bugfix for profiles output

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