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

Last change on this file since 4178 was 4131, checked in by monakurppa, 2 years ago

Several changes in the salsa aerosol module:

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