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

Last change on this file since 3721 was 3676, checked in by knoop, 6 years ago

Added comment with reference to PGI Compiler bug-ticket to the MERGE workaround in flow_statistics

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