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

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

Changed the name specific humidity to mixing ratio

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