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

Last change on this file since 3274 was 3274, checked in by knoop, 3 years ago

Modularization of all bulk cloud physics code components

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