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

Last change on this file since 2101 was 2101, checked in by suehring, 5 years ago

last commit documented

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