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

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

diagnostic output: Modularize diagnostic output, rename subroutines; formatting adjustments; allocate arrays only when required; add output of uu, vv, ww to enable variance calculation via temporal EC method; radiation: bugfix in masked data output; flow_statistics: Correct conversion to kinematic vertical scalar fluxes in case of pw-scheme and statistic regions

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