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

Last change on this file since 3662 was 3658, checked in by knoop, 5 years ago

OpenACC: flow_statistics partly ported to GPU

  • Property svn:keywords set to Id
File size: 111.6 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 3658 2019-01-07 20:28:54Z suehring $
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                   ! ki = MERGE( -1, 0, l == 0 )
1082                   ki = -1 + l
1083                   IF ( surf_def_h(l)%ns >= 1 )  THEN
1084                      DO  m = surf_def_h(l)%start_index(j,i),                  &
1085                              surf_def_h(l)%end_index(j,i)
1086                         k = surf_def_h(l)%k(m)
1087
1088                         !$ACC ATOMIC
1089                         sums_l(k+ki,12,tn) = sums_l(k+ki,12,tn) + &
1090                                    momentumflux_output_conversion(k+ki) * &
1091                                    surf_def_h(l)%usws(m) * rmask(j,i,sr)     ! w"u"
1092                         !$ACC ATOMIC
1093                         sums_l(k+ki,14,tn) = sums_l(k+ki,14,tn) + &
1094                                    momentumflux_output_conversion(k+ki) * &
1095                                    surf_def_h(l)%vsws(m) * rmask(j,i,sr)     ! w"v"
1096                         !$ACC ATOMIC
1097                         sums_l(k+ki,16,tn) = sums_l(k+ki,16,tn) + &
1098                                    heatflux_output_conversion(k+ki) * &
1099                                    surf_def_h(l)%shf(m)  * rmask(j,i,sr)     ! w"pt"
1100#if 0
1101                         sums_l(k+ki,58,tn) = sums_l(k+ki,58,tn) + &
1102                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
1103                         sums_l(k+ki,61,tn) = sums_l(k+ki,61,tn) + &
1104                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
1105#endif
1106#ifndef _OPENACC
1107                         IF ( ocean_mode )  THEN
1108                            sums_l(k+ki,65,tn) = sums_l(k+ki,65,tn) + &
1109                                       surf_def_h(l)%sasws(m) * rmask(j,i,sr)  ! w"sa"
1110                         ENDIF
1111                         IF ( humidity )  THEN
1112                            sums_l(k+ki,48,tn) = sums_l(k+ki,48,tn) +                     &
1113                                       waterflux_output_conversion(k+ki) *      &
1114                                       surf_def_h(l)%qsws(m) * rmask(j,i,sr)  ! w"q" (w"qv")
1115                            sums_l(k+ki,45,tn) = sums_l(k+ki,45,tn) + (                   &
1116                                       ( 1.0_wp + 0.61_wp * q(k+ki,j,i) ) *     &
1117                                       surf_def_h(l)%shf(m) + 0.61_wp * pt(k+ki,j,i) *      &
1118                                                  surf_def_h(l)%qsws(m) )                  &
1119                                       * heatflux_output_conversion(k+ki)
1120                            IF ( cloud_droplets )  THEN
1121                               sums_l(k+ki,45,tn) = sums_l(k+ki,45,tn) + (                &
1122                                         ( 1.0_wp + 0.61_wp * q(k+ki,j,i) -     &
1123                                           ql(k+ki,j,i) ) * surf_def_h(l)%shf(m) +          &
1124                                           0.61_wp * pt(k+ki,j,i) * surf_def_h(l)%qsws(m) ) &
1125                                          * heatflux_output_conversion(k+ki)
1126                            ENDIF
1127                            IF ( bulk_cloud_model )  THEN
1128!
1129!--                            Formula does not work if ql(k+ki) /= 0.0
1130                               sums_l(k+ki,51,tn) = sums_l(k+ki,51,tn) +                  &
1131                                          waterflux_output_conversion(k+ki) *   &
1132                                          surf_def_h(l)%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
1133                            ENDIF
1134                         ENDIF
1135                         IF ( passive_scalar )  THEN
1136                            sums_l(k+ki,117,tn) = sums_l(k+ki,117,tn) +                     &
1137                                        surf_def_h(l)%ssws(m) * rmask(j,i,sr) ! w"s"
1138                         ENDIF
1139#endif
1140
1141                      ENDDO
1142
1143                   ENDIF
1144                ENDDO
1145                IF ( surf_lsm_h%end_index(j,i) >=                              &
1146                     surf_lsm_h%start_index(j,i) )  THEN
1147                   m = surf_lsm_h%start_index(j,i)
1148                   !$ACC ATOMIC
1149                   sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
1150                                    momentumflux_output_conversion(nzb) * &
1151                                    surf_lsm_h%usws(m) * rmask(j,i,sr)     ! w"u"
1152                   !$ACC ATOMIC
1153                   sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
1154                                    momentumflux_output_conversion(nzb) * &
1155                                    surf_lsm_h%vsws(m) * rmask(j,i,sr)     ! w"v"
1156                   !$ACC ATOMIC
1157                   sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
1158                                    heatflux_output_conversion(nzb) * &
1159                                    surf_lsm_h%shf(m)  * rmask(j,i,sr)     ! w"pt"
1160#if 0
1161                   sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
1162                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
1163                   sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
1164                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
1165#endif
1166#ifndef _OPENACC
1167                   IF ( ocean_mode )  THEN
1168                      sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
1169                                       surf_lsm_h%sasws(m) * rmask(j,i,sr)  ! w"sa"
1170                   ENDIF
1171                   IF ( humidity )  THEN
1172                      sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
1173                                       waterflux_output_conversion(nzb) *      &
1174                                       surf_lsm_h%qsws(m) * rmask(j,i,sr)  ! w"q" (w"qv")
1175                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                   &
1176                                       ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) *     &
1177                                       surf_lsm_h%shf(m) + 0.61_wp * pt(nzb,j,i) *      &
1178                                                  surf_lsm_h%qsws(m) )                  &
1179                                       * heatflux_output_conversion(nzb)
1180                      IF ( cloud_droplets )  THEN
1181                         sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                &
1182                                         ( 1.0_wp + 0.61_wp * q(nzb,j,i) -     &
1183                                           ql(nzb,j,i) ) * surf_lsm_h%shf(m) +          &
1184                                           0.61_wp * pt(nzb,j,i) * surf_lsm_h%qsws(m) ) &
1185                                          * heatflux_output_conversion(nzb)
1186                      ENDIF
1187                      IF ( bulk_cloud_model )  THEN
1188!
1189!--                      Formula does not work if ql(nzb) /= 0.0
1190                         sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                  &
1191                                          waterflux_output_conversion(nzb) *   &
1192                                          surf_lsm_h%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
1193                      ENDIF
1194                   ENDIF
1195                   IF ( passive_scalar )  THEN
1196                      sums_l(nzb,117,tn) = sums_l(nzb,117,tn) +                     &
1197                                        surf_lsm_h%ssws(m) * rmask(j,i,sr) ! w"s"
1198                   ENDIF
1199#endif
1200
1201                ENDIF
1202                IF ( surf_usm_h%end_index(j,i) >=                              &
1203                     surf_usm_h%start_index(j,i) )  THEN
1204                   m = surf_usm_h%start_index(j,i)
1205                   !$ACC ATOMIC
1206                   sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
1207                                    momentumflux_output_conversion(nzb) * &
1208                                    surf_usm_h%usws(m) * rmask(j,i,sr)     ! w"u"
1209                   !$ACC ATOMIC
1210                   sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
1211                                    momentumflux_output_conversion(nzb) * &
1212                                    surf_usm_h%vsws(m) * rmask(j,i,sr)     ! w"v"
1213                   !$ACC ATOMIC
1214                   sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
1215                                    heatflux_output_conversion(nzb) * &
1216                                    surf_usm_h%shf(m)  * rmask(j,i,sr)     ! w"pt"
1217#if 0
1218                   sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
1219                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
1220                   sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
1221                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
1222#endif
1223#ifndef _OPENACC
1224                   IF ( ocean_mode )  THEN
1225                      sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
1226                                       surf_usm_h%sasws(m) * rmask(j,i,sr)  ! w"sa"
1227                   ENDIF
1228                   IF ( humidity )  THEN
1229                      sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
1230                                       waterflux_output_conversion(nzb) *      &
1231                                       surf_usm_h%qsws(m) * rmask(j,i,sr)  ! w"q" (w"qv")
1232                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                   &
1233                                       ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) *     &
1234                                       surf_usm_h%shf(m) + 0.61_wp * pt(nzb,j,i) *      &
1235                                                  surf_usm_h%qsws(m) )                  &
1236                                       * heatflux_output_conversion(nzb)
1237                      IF ( cloud_droplets )  THEN
1238                         sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                &
1239                                         ( 1.0_wp + 0.61_wp * q(nzb,j,i) -     &
1240                                           ql(nzb,j,i) ) * surf_usm_h%shf(m) +          &
1241                                           0.61_wp * pt(nzb,j,i) * surf_usm_h%qsws(m) ) &
1242                                          * heatflux_output_conversion(nzb)
1243                      ENDIF
1244                      IF ( bulk_cloud_model )  THEN
1245!
1246!--                      Formula does not work if ql(nzb) /= 0.0
1247                         sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                  &
1248                                          waterflux_output_conversion(nzb) *   &
1249                                          surf_usm_h%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
1250                      ENDIF
1251                   ENDIF
1252                   IF ( passive_scalar )  THEN
1253                      sums_l(nzb,117,tn) = sums_l(nzb,117,tn) +                     &
1254                                        surf_usm_h%ssws(m) * rmask(j,i,sr) ! w"s"
1255                   ENDIF
1256#endif
1257
1258                ENDIF
1259
1260             ENDIF
1261
1262#ifndef _OPENACC
1263             IF ( .NOT. neutral )  THEN
1264                IF ( surf_def_h(0)%end_index(j,i) >=                           &
1265                     surf_def_h(0)%start_index(j,i) )  THEN
1266                   m = surf_def_h(0)%start_index(j,i)
1267                   sums_l(nzb,112,tn) = sums_l(nzb,112,tn) +                   &
1268                                        surf_def_h(0)%ol(m)  * rmask(j,i,sr) ! L
1269                ENDIF
1270                IF ( surf_lsm_h%end_index(j,i) >=                              &
1271                     surf_lsm_h%start_index(j,i) )  THEN
1272                   m = surf_lsm_h%start_index(j,i)
1273                   sums_l(nzb,112,tn) = sums_l(nzb,112,tn) +                   &
1274                                        surf_lsm_h%ol(m)  * rmask(j,i,sr) ! L
1275                ENDIF
1276                IF ( surf_usm_h%end_index(j,i) >=                              &
1277                     surf_usm_h%start_index(j,i) )  THEN
1278                   m = surf_usm_h%start_index(j,i)
1279                   sums_l(nzb,112,tn) = sums_l(nzb,112,tn) +                   &
1280                                        surf_usm_h%ol(m)  * rmask(j,i,sr) ! L
1281                ENDIF
1282             ENDIF
1283
1284             IF ( radiation )  THEN
1285                IF ( surf_def_h(0)%end_index(j,i) >=                           &
1286                     surf_def_h(0)%start_index(j,i) )  THEN
1287                   m = surf_def_h(0)%start_index(j,i)
1288                   sums_l(nzb,99,tn)  = sums_l(nzb,99,tn)  +                   &
1289                                        surf_def_h(0)%rad_net(m) * rmask(j,i,sr)
1290                   sums_l(nzb,100,tn) = sums_l(nzb,100,tn)  +                  &
1291                                        surf_def_h(0)%rad_lw_in(m) * rmask(j,i,sr)
1292                   sums_l(nzb,101,tn) = sums_l(nzb,101,tn)  +                  &
1293                                        surf_def_h(0)%rad_lw_out(m) * rmask(j,i,sr)
1294                   sums_l(nzb,102,tn) = sums_l(nzb,102,tn)  +                  &
1295                                        surf_def_h(0)%rad_sw_in(m) * rmask(j,i,sr)
1296                   sums_l(nzb,103,tn) = sums_l(nzb,103,tn)  +                  &
1297                                        surf_def_h(0)%rad_sw_out(m) * rmask(j,i,sr)
1298                ENDIF
1299                IF ( surf_lsm_h%end_index(j,i) >=                              &
1300                     surf_lsm_h%start_index(j,i) )  THEN
1301                   m = surf_lsm_h%start_index(j,i)
1302                   sums_l(nzb,99,tn)  = sums_l(nzb,99,tn)  +                   &
1303                                        surf_lsm_h%rad_net(m) * rmask(j,i,sr)
1304                   sums_l(nzb,100,tn) = sums_l(nzb,100,tn)  +                  &
1305                                        surf_lsm_h%rad_lw_in(m) * rmask(j,i,sr)
1306                   sums_l(nzb,101,tn) = sums_l(nzb,101,tn)  +                  &
1307                                        surf_lsm_h%rad_lw_out(m) * rmask(j,i,sr)
1308                   sums_l(nzb,102,tn) = sums_l(nzb,102,tn)  +                  &
1309                                        surf_lsm_h%rad_sw_in(m) * rmask(j,i,sr)
1310                   sums_l(nzb,103,tn) = sums_l(nzb,103,tn)  +                  &
1311                                        surf_lsm_h%rad_sw_out(m) * rmask(j,i,sr)
1312                ENDIF
1313                IF ( surf_usm_h%end_index(j,i) >=                              &
1314                     surf_usm_h%start_index(j,i) )  THEN
1315                   m = surf_usm_h%start_index(j,i)
1316                   sums_l(nzb,99,tn)  = sums_l(nzb,99,tn)  +                   &
1317                                        surf_usm_h%rad_net(m) * rmask(j,i,sr)
1318                   sums_l(nzb,100,tn) = sums_l(nzb,100,tn)  +                  &
1319                                        surf_usm_h%rad_lw_in(m) * rmask(j,i,sr)
1320                   sums_l(nzb,101,tn) = sums_l(nzb,101,tn)  +                  &
1321                                        surf_usm_h%rad_lw_out(m) * rmask(j,i,sr)
1322                   sums_l(nzb,102,tn) = sums_l(nzb,102,tn)  +                  &
1323                                        surf_usm_h%rad_sw_in(m) * rmask(j,i,sr)
1324                   sums_l(nzb,103,tn) = sums_l(nzb,103,tn)  +                  &
1325                                        surf_usm_h%rad_sw_out(m) * rmask(j,i,sr)
1326                ENDIF
1327
1328#if defined ( __rrtmg )
1329                IF ( radiation_scheme == 'rrtmg' )  THEN
1330
1331                   IF ( surf_def_h(0)%end_index(j,i) >=                        &
1332                        surf_def_h(0)%start_index(j,i) )  THEN
1333                      m = surf_def_h(0)%start_index(j,i)
1334                      sums_l(nzb,108,tn)  = sums_l(nzb,108,tn)  +              &
1335                                   surf_def_h(0)%rrtm_aldif(0,m) * rmask(j,i,sr)
1336                      sums_l(nzb,109,tn) = sums_l(nzb,109,tn)  +               &
1337                                   surf_def_h(0)%rrtm_aldir(0,m) * rmask(j,i,sr)
1338                      sums_l(nzb,110,tn) = sums_l(nzb,110,tn)  +               &
1339                                   surf_def_h(0)%rrtm_asdif(0,m) * rmask(j,i,sr)
1340                      sums_l(nzb,111,tn) = sums_l(nzb,111,tn)  +               &
1341                                   surf_def_h(0)%rrtm_asdir(0,m) * rmask(j,i,sr)
1342                   ENDIF
1343                   IF ( surf_lsm_h%end_index(j,i) >=                           &
1344                        surf_lsm_h%start_index(j,i) )  THEN
1345                      m = surf_lsm_h%start_index(j,i)
1346                      sums_l(nzb,108,tn)  = sums_l(nzb,108,tn)  +              &
1347                               SUM( surf_lsm_h%frac(:,m) *                     &
1348                                    surf_lsm_h%rrtm_aldif(:,m) ) * rmask(j,i,sr)
1349                      sums_l(nzb,109,tn) = sums_l(nzb,109,tn)  +               &
1350                               SUM( surf_lsm_h%frac(:,m) *                     &
1351                                    surf_lsm_h%rrtm_aldir(:,m) ) * rmask(j,i,sr)
1352                      sums_l(nzb,110,tn) = sums_l(nzb,110,tn)  +               &
1353                               SUM( surf_lsm_h%frac(:,m) *                     &
1354                                    surf_lsm_h%rrtm_asdif(:,m) ) * rmask(j,i,sr)
1355                      sums_l(nzb,111,tn) = sums_l(nzb,111,tn)  +               &
1356                               SUM( surf_lsm_h%frac(:,m) *                     &
1357                                    surf_lsm_h%rrtm_asdir(:,m) ) * rmask(j,i,sr)
1358                   ENDIF
1359                   IF ( surf_usm_h%end_index(j,i) >=                           &
1360                        surf_usm_h%start_index(j,i) )  THEN
1361                      m = surf_usm_h%start_index(j,i)
1362                      sums_l(nzb,108,tn)  = sums_l(nzb,108,tn)  +              &
1363                               SUM( surf_usm_h%frac(:,m) *                     &
1364                                    surf_usm_h%rrtm_aldif(:,m) ) * rmask(j,i,sr)
1365                      sums_l(nzb,109,tn) = sums_l(nzb,109,tn)  +               &
1366                               SUM( surf_usm_h%frac(:,m) *                     &
1367                                    surf_usm_h%rrtm_aldir(:,m) ) * rmask(j,i,sr)
1368                      sums_l(nzb,110,tn) = sums_l(nzb,110,tn)  +               &
1369                               SUM( surf_usm_h%frac(:,m) *                     &
1370                                    surf_usm_h%rrtm_asdif(:,m) ) * rmask(j,i,sr)
1371                      sums_l(nzb,111,tn) = sums_l(nzb,111,tn)  +               &
1372                               SUM( surf_usm_h%frac(:,m) *                     &
1373                                    surf_usm_h%rrtm_asdir(:,m) ) * rmask(j,i,sr)
1374                   ENDIF
1375
1376                ENDIF
1377#endif
1378             ENDIF
1379#endif
1380!
1381!--          Subgridscale fluxes at the top surface
1382             IF ( use_top_fluxes )  THEN
1383                m = surf_def_h(2)%start_index(j,i)
1384                !$ACC ATOMIC
1385                sums_l(nzt,12,tn) = sums_l(nzt,12,tn) + &
1386                                    momentumflux_output_conversion(nzt) * &
1387                                    surf_def_h(2)%usws(m) * rmask(j,i,sr)    ! w"u"
1388                !$ACC ATOMIC
1389                sums_l(nzt+1,12,tn) = sums_l(nzt+1,12,tn) + &
1390                                    momentumflux_output_conversion(nzt+1) * &
1391                                    surf_def_h(2)%usws(m) * rmask(j,i,sr)    ! w"u"
1392                !$ACC ATOMIC
1393                sums_l(nzt,14,tn) = sums_l(nzt,14,tn) + &
1394                                    momentumflux_output_conversion(nzt) * &
1395                                    surf_def_h(2)%vsws(m) * rmask(j,i,sr)    ! w"v"
1396                !$ACC ATOMIC
1397                sums_l(nzt+1,14,tn) = sums_l(nzt+1,14,tn) + &
1398                                    momentumflux_output_conversion(nzt+1) * &
1399                                    surf_def_h(2)%vsws(m) * rmask(j,i,sr)    ! w"v"
1400                !$ACC ATOMIC
1401                sums_l(nzt,16,tn) = sums_l(nzt,16,tn) + &
1402                                    heatflux_output_conversion(nzt) * &
1403                                    surf_def_h(2)%shf(m)  * rmask(j,i,sr)   ! w"pt"
1404                !$ACC ATOMIC
1405                sums_l(nzt+1,16,tn) = sums_l(nzt+1,16,tn) + &
1406                                    heatflux_output_conversion(nzt+1) * &
1407                                    surf_def_h(2)%shf(m)  * rmask(j,i,sr)   ! w"pt"
1408#if 0
1409                sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + &
1410                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
1411                sums_l(nzt:nzt+1,61,tn) = sums_l(nzt:nzt+1,61,tn) + &
1412                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
1413#endif
1414#ifndef _OPENACC
1415                IF ( ocean_mode )  THEN
1416                   sums_l(nzt,65,tn) = sums_l(nzt,65,tn) + &
1417                                       surf_def_h(2)%sasws(m) * rmask(j,i,sr)  ! w"sa"
1418                ENDIF
1419                IF ( humidity )  THEN
1420                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) +                     &
1421                                       waterflux_output_conversion(nzt) *      &
1422                                       surf_def_h(2)%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
1423                   sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                   &
1424                                       ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) *     &
1425                                       surf_def_h(2)%shf(m) +                  &
1426                                       0.61_wp * pt(nzt,j,i) *    &
1427                                       surf_def_h(2)%qsws(m) )      &
1428                                       * heatflux_output_conversion(nzt)
1429                   IF ( cloud_droplets )  THEN
1430                      sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                &
1431                                          ( 1.0_wp + 0.61_wp * q(nzt,j,i) -    &
1432                                            ql(nzt,j,i) ) *                    &
1433                                            surf_def_h(2)%shf(m) +             &
1434                                           0.61_wp * pt(nzt,j,i) *             &
1435                                           surf_def_h(2)%qsws(m) )&
1436                                           * heatflux_output_conversion(nzt)
1437                   ENDIF
1438                   IF ( bulk_cloud_model )  THEN
1439!
1440!--                   Formula does not work if ql(nzb) /= 0.0
1441                      sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + &   ! w"q" (w"qv")
1442                                          waterflux_output_conversion(nzt) *   &
1443                                          surf_def_h(2)%qsws(m) * rmask(j,i,sr)
1444                   ENDIF
1445                ENDIF
1446                IF ( passive_scalar )  THEN
1447                   sums_l(nzt,117,tn) = sums_l(nzt,117,tn) + &
1448                                        surf_def_h(2)%ssws(m) * rmask(j,i,sr) ! w"s"
1449                ENDIF
1450#endif
1451             ENDIF
1452
1453!
1454!--          Resolved fluxes (can be computed for all horizontal points)
1455!--          NOTE: for simplicity, nzb_s_inner is used below, although strictly
1456!--          ----  speaking the following k-loop would have to be split up and
1457!--                rearranged according to the staggered grid.
1458             DO  k = nzb, nzt
1459                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
1460                ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +                  &
1461                                 u(k+1,j,i) - hom(k+1,1,1,sr) )
1462                vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +                  &
1463                                 v(k+1,j,i) - hom(k+1,1,2,sr) )
1464                pts = 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) +                 &
1465                                 pt(k+1,j,i) - hom(k+1,1,4,sr) )
1466
1467!--             Higher moments
1468                !$ACC ATOMIC
1469                sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 *        &
1470                                                    rmask(j,i,sr) * flag
1471                !$ACC ATOMIC
1472                sums_l(k,36,tn) = sums_l(k,36,tn) + pts**2 * w(k,j,i) *        &
1473                                                    rmask(j,i,sr) * flag
1474
1475!
1476!--             Salinity flux and density (density does not belong to here,
1477!--             but so far there is no other suitable place to calculate)
1478#ifndef _OPENACC
1479                IF ( ocean_mode )  THEN
1480                   IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
1481                      pts = 0.5_wp * ( sa(k,j,i)   - hom(k,1,23,sr) +          &
1482                                       sa(k+1,j,i) - hom(k+1,1,23,sr) )
1483                      sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) *     &
1484                                        rmask(j,i,sr) * flag
1485                   ENDIF
1486                   sums_l(k,64,tn) = sums_l(k,64,tn) + rho_ocean(k,j,i) *      &
1487                                                       rmask(j,i,sr) * flag
1488                   sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) *           &
1489                                                       rmask(j,i,sr) * flag
1490                ENDIF
1491
1492!
1493!--             Buoyancy flux, water flux, humidity flux, liquid water
1494!--             content, rain drop concentration and rain water content
1495                IF ( humidity )  THEN
1496                   IF ( bulk_cloud_model .OR. cloud_droplets )  THEN
1497                      pts = 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +         &
1498                                    vpt(k+1,j,i) - hom(k+1,1,44,sr) )
1499                      sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *     &
1500                                               heatflux_output_conversion(k) * &
1501                                                          rmask(j,i,sr) * flag
1502                      sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * rmask(j,i,sr) &
1503                                                                    * flag
1504
1505                      IF ( .NOT. cloud_droplets )  THEN
1506                         pts = 0.5_wp *                                        &
1507                              ( ( q(k,j,i) - ql(k,j,i) ) -                     &
1508                              hom(k,1,42,sr) +                                 &
1509                              ( q(k+1,j,i) - ql(k+1,j,i) ) -                   &
1510                              hom(k+1,1,42,sr) )
1511                         sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) *  &
1512                                             waterflux_output_conversion(k) *  &
1513                                                             rmask(j,i,sr)  *  &
1514                                                             flag
1515                         sums_l(k,75,tn) = sums_l(k,75,tn) + qc(k,j,i) *       &
1516                                                             rmask(j,i,sr) *   &
1517                                                             flag
1518                         sums_l(k,76,tn) = sums_l(k,76,tn) + prr(k,j,i) *      &
1519                                                             rmask(j,i,sr) *   &
1520                                                             flag
1521                         IF ( microphysics_morrison )  THEN
1522                            sums_l(k,123,tn) = sums_l(k,123,tn) + nc(k,j,i) *  &
1523                                                                rmask(j,i,sr) *&
1524                                                                flag
1525                         ENDIF
1526                         IF ( microphysics_seifert )  THEN
1527                            sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) *    &
1528                                                                rmask(j,i,sr) *&
1529                                                                flag 
1530                            sums_l(k,74,tn) = sums_l(k,74,tn) + qr(k,j,i) *    &
1531                                                                rmask(j,i,sr) *&
1532                                                                flag
1533                         ENDIF
1534                      ENDIF
1535
1536                   ELSE
1537                      IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
1538                         pts = 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +      &
1539                                          vpt(k+1,j,i) - hom(k+1,1,44,sr) )
1540                         sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *  &
1541                                              heatflux_output_conversion(k) *  &
1542                                                             rmask(j,i,sr)  *  &
1543                                                             flag
1544                      ELSE IF ( ws_scheme_sca .AND. sr == 0 )  THEN
1545                         sums_l(k,46,tn) = ( ( 1.0_wp + 0.61_wp *              & 
1546                                               hom(k,1,41,sr) ) *              &
1547                                             sums_l(k,17,tn) +                 &
1548                                             0.61_wp * hom(k,1,4,sr) *         &
1549                                             sums_l(k,49,tn)                   &
1550                                           ) * heatflux_output_conversion(k) * &
1551                                               flag
1552                      END IF
1553                   END IF
1554                ENDIF
1555!
1556!--             Passive scalar flux
1557                IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca                &
1558                     .OR. sr /= 0 ) )  THEN
1559                   pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,115,sr) +             &
1560                                    s(k+1,j,i) - hom(k+1,1,115,sr) )
1561                   sums_l(k,114,tn) = sums_l(k,114,tn) + pts * w(k,j,i) *      &
1562                                                       rmask(j,i,sr) * flag
1563                ENDIF
1564#endif
1565
1566!
1567!--             Energy flux w*e*
1568!--             has to be adjusted
1569                !$ACC ATOMIC
1570                sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5_wp *        &
1571                                             ( ust**2 + vst**2 + w(k,j,i)**2 ) &
1572                                           * rho_air_zw(k)                     &
1573                                           * momentumflux_output_conversion(k) &
1574                                           * rmask(j,i,sr) * flag
1575             ENDDO
1576          ENDDO
1577       ENDDO
1578       !$OMP END PARALLEL
1579
1580       !$ACC UPDATE &
1581       !$ACC HOST(sums_l(:,12,tn), sums_l(:,14,tn), sums_l(:,16,tn)) &
1582       !$ACC HOST(sums_l(:,35,tn), sums_l(:,36,tn), sums_l(:,37,tn))
1583!
1584!--    Treat land-surface quantities according to new wall model structure.
1585       IF ( land_surface )  THEN
1586          tn = 0
1587          !$OMP PARALLEL PRIVATE( i, j, m, tn )
1588          !$ tn = omp_get_thread_num()
1589          !$OMP DO
1590          DO  m = 1, surf_lsm_h%ns
1591             i = surf_lsm_h%i(m)
1592             j = surf_lsm_h%j(m)
1593       
1594             IF ( i >= nxl  .AND.  i <= nxr  .AND.                             &
1595                  j >= nys  .AND.  j <= nyn )  THEN
1596                sums_l(nzb,93,tn)  = sums_l(nzb,93,tn) + surf_lsm_h%ghf(m)
1597                sums_l(nzb,94,tn)  = sums_l(nzb,94,tn) + surf_lsm_h%qsws_liq(m)
1598                sums_l(nzb,95,tn)  = sums_l(nzb,95,tn) + surf_lsm_h%qsws_soil(m)
1599                sums_l(nzb,96,tn)  = sums_l(nzb,96,tn) + surf_lsm_h%qsws_veg(m)
1600                sums_l(nzb,97,tn)  = sums_l(nzb,97,tn) + surf_lsm_h%r_a(m)
1601                sums_l(nzb,98,tn) = sums_l(nzb,98,tn)+ surf_lsm_h%r_s(m)
1602             ENDIF
1603          ENDDO
1604          !$OMP END PARALLEL
1605
1606          tn = 0
1607          !$OMP PARALLEL PRIVATE( i, j, k, m, tn )
1608          !$ tn = omp_get_thread_num()
1609          !$OMP DO
1610          DO  m = 1, surf_lsm_h%ns
1611
1612             i = surf_lsm_h%i(m)           
1613             j = surf_lsm_h%j(m)
1614
1615             IF ( i >= nxl  .AND.  i <= nxr  .AND.                             &
1616                  j >= nys  .AND.  j <= nyn )  THEN
1617
1618                DO  k = nzb_soil, nzt_soil
1619                   sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil_h%var_2d(k,m)  &
1620                                      * rmask(j,i,sr)
1621                   sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil_h%var_2d(k,m)  &
1622                                      * rmask(j,i,sr)
1623                ENDDO
1624             ENDIF
1625          ENDDO
1626          !$OMP END PARALLEL
1627       ENDIF
1628!
1629!--    For speed optimization fluxes which have been computed in part directly
1630!--    inside the WS advection routines are treated seperatly
1631!--    Momentum fluxes first:
1632
1633       tn = 0
1634       !$OMP PARALLEL PRIVATE( i, j, k, tn, flag )
1635       !$ tn = omp_get_thread_num()
1636       IF ( .NOT. ws_scheme_mom .OR. sr /= 0  )  THEN
1637          !$OMP DO
1638          DO  i = nxl, nxr
1639             DO  j = nys, nyn
1640                DO  k = nzb, nzt
1641!
1642!--                Flag 23 is used to mask surface fluxes as well as model-top
1643!--                fluxes, which are added further below.
1644                   flag = MERGE( 1.0_wp, 0.0_wp,                               &
1645                                 BTEST( wall_flags_0(k,j,i), 23 ) ) *          &
1646                          MERGE( 1.0_wp, 0.0_wp,                               &
1647                                 BTEST( wall_flags_0(k,j,i), 9  ) )
1648
1649                   ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +               &
1650                                    u(k+1,j,i) - hom(k+1,1,1,sr) )
1651                   vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +               &
1652                                    v(k+1,j,i) - hom(k+1,1,2,sr) )
1653!
1654!--                Momentum flux w*u*
1655                   sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5_wp *                &
1656                                                     ( w(k,j,i-1) + w(k,j,i) ) &
1657                                           * rho_air_zw(k)                     &
1658                                           * momentumflux_output_conversion(k) &
1659                                                     * ust * rmask(j,i,sr)     &
1660                                                           * flag
1661!
1662!--                Momentum flux w*v*
1663                   sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5_wp *                &
1664                                                     ( w(k,j-1,i) + w(k,j,i) ) &
1665                                           * rho_air_zw(k)                     &
1666                                           * momentumflux_output_conversion(k) &
1667                                                     * vst * rmask(j,i,sr)     &
1668                                                           * flag
1669                ENDDO
1670             ENDDO
1671          ENDDO
1672
1673       ENDIF
1674       IF ( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
1675          !$OMP DO
1676          DO  i = nxl, nxr
1677             DO  j = nys, nyn
1678                DO  k = nzb, nzt
1679                   flag = MERGE( 1.0_wp, 0.0_wp,                               &
1680                                 BTEST( wall_flags_0(k,j,i), 23 ) ) *          &
1681                          MERGE( 1.0_wp, 0.0_wp,                               &
1682                                 BTEST( wall_flags_0(k,j,i), 9  ) )
1683!
1684!--                Vertical heat flux
1685                   sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5_wp *                &
1686                           ( pt(k,j,i)   - hom(k,1,4,sr) +                     &
1687                             pt(k+1,j,i) - hom(k+1,1,4,sr) )                   &
1688                           * heatflux_output_conversion(k)                     &
1689                           * w(k,j,i) * rmask(j,i,sr) * flag
1690                   IF ( humidity )  THEN
1691                      pts = 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +           &
1692                                      q(k+1,j,i) - hom(k+1,1,41,sr) )
1693                      sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *     &
1694                                       waterflux_output_conversion(k) *        &
1695                                       rmask(j,i,sr) * flag
1696                   ENDIF
1697                   IF ( passive_scalar )  THEN
1698                      pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,115,sr) +           &
1699                                      s(k+1,j,i) - hom(k+1,1,115,sr) )
1700                      sums_l(k,114,tn) = sums_l(k,114,tn) + pts * w(k,j,i) *    &
1701                                        rmask(j,i,sr) * flag
1702                   ENDIF
1703                ENDDO
1704             ENDDO
1705          ENDDO
1706
1707       ENDIF
1708
1709!
1710!--    Density at top follows Neumann condition
1711       IF ( ocean_mode )  THEN
1712          sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
1713          sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
1714       ENDIF
1715
1716!
1717!--    Divergence of vertical flux of resolved scale energy and pressure
1718!--    fluctuations as well as flux of pressure fluctuation itself (68).
1719!--    First calculate the products, then the divergence.
1720!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
1721       IF ( hom(nzb+1,2,55,0) /= 0.0_wp  .OR.  hom(nzb+1,2,68,0) /= 0.0_wp )   &
1722       THEN
1723          sums_ll = 0.0_wp  ! local array
1724
1725          !$OMP DO
1726          DO  i = nxl, nxr
1727             DO  j = nys, nyn
1728                DO  k = nzb+1, nzt
1729                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1730
1731                   sums_ll(k,1) = sums_ll(k,1) + 0.5_wp * w(k,j,i) * (         &
1732                  ( 0.25_wp * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1) )  &
1733                            - 0.5_wp * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) ) )**2&
1734                + ( 0.25_wp * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i) )  &
1735                            - 0.5_wp * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) ) )**2&
1736                + w(k,j,i)**2                                        ) * flag
1737
1738                   sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i)             &
1739                                       * ( ( p(k,j,i) + p(k+1,j,i) )           &
1740                                         / momentumflux_output_conversion(k) ) &
1741                                       * flag
1742
1743                ENDDO
1744             ENDDO
1745          ENDDO
1746          sums_ll(0,1)     = 0.0_wp    ! because w is zero at the bottom
1747          sums_ll(nzt+1,1) = 0.0_wp
1748          sums_ll(0,2)     = 0.0_wp
1749          sums_ll(nzt+1,2) = 0.0_wp
1750
1751          DO  k = nzb+1, nzt
1752             sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
1753             sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
1754             sums_l(k,68,tn) = sums_ll(k,2)
1755          ENDDO
1756          sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
1757          sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
1758          sums_l(nzb,68,tn) = 0.0_wp    ! because w* = 0 at nzb
1759
1760       ENDIF
1761
1762!
1763!--    Divergence of vertical flux of SGS TKE and the flux itself (69)
1764       IF ( hom(nzb+1,2,57,0) /= 0.0_wp  .OR.  hom(nzb+1,2,69,0) /= 0.0_wp )   &
1765       THEN
1766          !$OMP DO
1767          DO  i = nxl, nxr
1768             DO  j = nys, nyn
1769                DO  k = nzb+1, nzt
1770
1771                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1772
1773                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5_wp * (              &
1774                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
1775                 - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
1776                                                                ) * ddzw(k)    &
1777                                                                  * flag
1778
1779                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5_wp * (              &
1780                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
1781                                                                )  * flag
1782
1783                ENDDO
1784             ENDDO
1785          ENDDO
1786          sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
1787          sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
1788
1789       ENDIF
1790
1791!
1792!--    Horizontal heat fluxes (subgrid, resolved, total).
1793!--    Do it only, if profiles shall be plotted.
1794       IF ( hom(nzb+1,2,58,0) /= 0.0_wp ) THEN
1795
1796          !$OMP DO
1797          DO  i = nxl, nxr
1798             DO  j = nys, nyn
1799                DO  k = nzb+1, nzt
1800                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1801!
1802!--                Subgrid horizontal heat fluxes u"pt", v"pt"
1803                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5_wp *                &
1804                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
1805                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
1806                                               * rho_air_zw(k)                 &
1807                                               * heatflux_output_conversion(k) &
1808                                                 * ddx * rmask(j,i,sr) * flag
1809                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
1810                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
1811                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
1812                                               * rho_air_zw(k)                 &
1813                                               * heatflux_output_conversion(k) &
1814                                                 * ddy * rmask(j,i,sr) * flag
1815!
1816!--                Resolved horizontal heat fluxes u*pt*, v*pt*
1817                   sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
1818                                                  ( u(k,j,i) - hom(k,1,1,sr) ) &
1819                                    * 0.5_wp * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
1820                                                 pt(k,j,i)   - hom(k,1,4,sr) ) &
1821                                               * heatflux_output_conversion(k) &
1822                                               * flag
1823                   pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
1824                                    pt(k,j,i)   - hom(k,1,4,sr) )
1825                   sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
1826                                                  ( v(k,j,i) - hom(k,1,2,sr) ) &
1827                                    * 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
1828                                                 pt(k,j,i)   - hom(k,1,4,sr) ) &
1829                                               * heatflux_output_conversion(k) &
1830                                               * flag
1831                ENDDO
1832             ENDDO
1833          ENDDO
1834!
1835!--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
1836          sums_l(nzb,58,tn) = 0.0_wp
1837          sums_l(nzb,59,tn) = 0.0_wp
1838          sums_l(nzb,60,tn) = 0.0_wp
1839          sums_l(nzb,61,tn) = 0.0_wp
1840          sums_l(nzb,62,tn) = 0.0_wp
1841          sums_l(nzb,63,tn) = 0.0_wp
1842
1843       ENDIF
1844       !$OMP END PARALLEL
1845
1846!
1847!--    Collect current large scale advection and subsidence tendencies for
1848!--    data output
1849       IF ( large_scale_forcing  .AND.  ( simulated_time > 0.0_wp ) )  THEN
1850!
1851!--       Interpolation in time of LSF_DATA
1852          nt = 1
1853          DO WHILE ( simulated_time - dt_3d > time_vert(nt) )
1854             nt = nt + 1
1855          ENDDO
1856          IF ( simulated_time - dt_3d /= time_vert(nt) )  THEN
1857            nt = nt - 1
1858          ENDIF
1859
1860          fac = ( simulated_time - dt_3d - time_vert(nt) )                     &
1861                / ( time_vert(nt+1)-time_vert(nt) )
1862
1863
1864          DO  k = nzb, nzt
1865             sums_ls_l(k,0) = td_lsa_lpt(k,nt)                                 &
1866                              + fac * ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )
1867             sums_ls_l(k,1) = td_lsa_q(k,nt)                                   &
1868                              + fac * ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) )
1869          ENDDO
1870
1871          sums_ls_l(nzt+1,0) = sums_ls_l(nzt,0)
1872          sums_ls_l(nzt+1,1) = sums_ls_l(nzt,1)
1873
1874          IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
1875
1876             DO  k = nzb, nzt
1877                sums_ls_l(k,2) = td_sub_lpt(k,nt) + fac *                      &
1878                                 ( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )
1879                sums_ls_l(k,3) = td_sub_q(k,nt) + fac *                        &
1880                                 ( td_sub_q(k,nt+1) - td_sub_q(k,nt) )
1881             ENDDO
1882
1883             sums_ls_l(nzt+1,2) = sums_ls_l(nzt,2)
1884             sums_ls_l(nzt+1,3) = sums_ls_l(nzt,3)
1885
1886          ENDIF
1887
1888       ENDIF
1889
1890       tn = 0
1891       !$OMP PARALLEL PRIVATE( i, j, k, tn )
1892       !$ tn = omp_get_thread_num()       
1893       IF ( radiation .AND. radiation_scheme == 'rrtmg' )  THEN
1894          !$OMP DO
1895          DO  i = nxl, nxr
1896             DO  j =  nys, nyn
1897                DO  k = nzb+1, nzt+1
1898                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1899
1900                   sums_l(k,100,tn)  = sums_l(k,100,tn)  + rad_lw_in(k,j,i)    &
1901                                       * rmask(j,i,sr) * flag
1902                   sums_l(k,101,tn)  = sums_l(k,101,tn)  + rad_lw_out(k,j,i)   &
1903                                       * rmask(j,i,sr) * flag
1904                   sums_l(k,102,tn)  = sums_l(k,102,tn)  + rad_sw_in(k,j,i)    &
1905                                       * rmask(j,i,sr) * flag
1906                   sums_l(k,103,tn)  = sums_l(k,103,tn)  + rad_sw_out(k,j,i)   &
1907                                       * rmask(j,i,sr) * flag
1908                   sums_l(k,104,tn)  = sums_l(k,104,tn)  + rad_lw_cs_hr(k,j,i) &
1909                                       * rmask(j,i,sr) * flag
1910                   sums_l(k,105,tn)  = sums_l(k,105,tn)  + rad_lw_hr(k,j,i)    &
1911                                       * rmask(j,i,sr) * flag
1912                   sums_l(k,106,tn)  = sums_l(k,106,tn)  + rad_sw_cs_hr(k,j,i) &
1913                                       * rmask(j,i,sr) * flag
1914                   sums_l(k,107,tn)  = sums_l(k,107,tn)  + rad_sw_hr(k,j,i)    &
1915                                       * rmask(j,i,sr) * flag
1916                ENDDO
1917             ENDDO
1918          ENDDO
1919       ENDIF
1920
1921!
1922!--    Calculate the profiles for all other modules
1923       CALL module_interface_statistics( 'profiles', sr, tn, dots_max )
1924       !$OMP END PARALLEL
1925
1926!
1927!--    Summation of thread sums
1928       IF ( threads_per_task > 1 )  THEN
1929          DO  i = 1, threads_per_task-1
1930             sums_l(:,3,0)          = sums_l(:,3,0) + sums_l(:,3,i)
1931             sums_l(:,4:40,0)       = sums_l(:,4:40,0) + sums_l(:,4:40,i)
1932             sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
1933                                      sums_l(:,45:pr_palm,i)
1934             IF ( max_pr_user > 0 )  THEN
1935                sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
1936                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
1937                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
1938             ENDIF
1939
1940             IF ( air_chemistry )  THEN
1941                IF ( max_pr_cs > 0 )  THEN                                 
1942                     sums_l(:,pr_palm+max_pr_user+1:pr_palm + max_pr_user+ max_pr_cs,0) =          &
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,i)
1945
1946                ENDIF
1947             ENDIF
1948          ENDDO
1949       ENDIF
1950
1951#if defined( __parallel )
1952
1953!
1954!--    Compute total sum from local sums
1955       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1956       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL,   &
1957                           MPI_SUM, comm2d, ierr )
1958       IF ( large_scale_forcing )  THEN
1959          CALL MPI_ALLREDUCE( sums_ls_l(nzb,2), sums(nzb,83), ngp_sums_ls,     &
1960                              MPI_REAL, MPI_SUM, comm2d, ierr )
1961       ENDIF
1962
1963       IF ( air_chemistry  .AND.  max_pr_cs > 0 )  THEN
1964          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1965             DO  i = 1, max_pr_cs
1966                CALL MPI_ALLREDUCE( sums_l(nzb,pr_palm+max_pr_user+i,0),       &
1967                                    sums(nzb,pr_palm+max_pr_user+i),           &
1968                                    nzt+2-nzb, MPI_REAL, MPI_SUM, comm2d, ierr )
1969             ENDDO
1970       ENDIF
1971
1972#else
1973       sums = sums_l(:,:,0)
1974       IF ( large_scale_forcing )  THEN
1975          sums(:,81:88) = sums_ls_l
1976       ENDIF
1977#endif
1978
1979!
1980!--    Final values are obtained by division by the total number of grid points
1981!--    used for summation. After that store profiles.
1982!--    Check, if statistical regions do contain at least one grid point at the
1983!--    respective k-level, otherwise division by zero will lead to undefined
1984!--    values, which may cause e.g. problems with NetCDF output
1985!--    Profiles:
1986       DO  k = nzb, nzt+1
1987          sums(k,3)             = sums(k,3)             / ngp_2dh(sr)
1988          sums(k,12:22)         = sums(k,12:22)         / ngp_2dh(sr)
1989          sums(k,30:32)         = sums(k,30:32)         / ngp_2dh(sr)
1990          sums(k,35:39)         = sums(k,35:39)         / ngp_2dh(sr)
1991          sums(k,45:53)         = sums(k,45:53)         / ngp_2dh(sr)
1992          sums(k,55:63)         = sums(k,55:63)         / ngp_2dh(sr)
1993          sums(k,81:88)         = sums(k,81:88)         / ngp_2dh(sr)
1994          sums(k,89:112)        = sums(k,89:112)        / ngp_2dh(sr)
1995          sums(k,114)           = sums(k,114)           / ngp_2dh(sr)
1996          sums(k,117)           = sums(k,117)           / ngp_2dh(sr)
1997          IF ( ngp_2dh_s_inner(k,sr) /= 0 )  THEN
1998             sums(k,8:11)          = sums(k,8:11)          / ngp_2dh_s_inner(k,sr)
1999             sums(k,23:29)         = sums(k,23:29)         / ngp_2dh_s_inner(k,sr)
2000             sums(k,33:34)         = sums(k,33:34)         / ngp_2dh_s_inner(k,sr)
2001             sums(k,40)            = sums(k,40)            / ngp_2dh_s_inner(k,sr)
2002             sums(k,54)            = sums(k,54)            / ngp_2dh_s_inner(k,sr)
2003             sums(k,64)            = sums(k,64)            / ngp_2dh_s_inner(k,sr)
2004             sums(k,70:80)         = sums(k,70:80)         / ngp_2dh_s_inner(k,sr)
2005             sums(k,116)           = sums(k,116)           / ngp_2dh_s_inner(k,sr)
2006             sums(k,118:pr_palm-2) = sums(k,118:pr_palm-2) / ngp_2dh_s_inner(k,sr)
2007             sums(k,123)           = sums(k,123) * ngp_2dh_s_inner(k,sr)  / ngp_2dh(sr)
2008          ENDIF
2009       ENDDO
2010
2011!--    u* and so on
2012!--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
2013!--    size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
2014!--    above the topography, they are being divided by ngp_2dh(sr)
2015       sums(nzb:nzb+3,pr_palm)    = sums(nzb:nzb+3,pr_palm)    / &
2016                                    ngp_2dh(sr)
2017       sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)       / &    ! qs
2018                                    ngp_2dh(sr)
2019       sums(nzb+13,pr_palm)       = sums(nzb+13,pr_palm)       / &    ! ss
2020                                    ngp_2dh(sr)
2021       sums(nzb+14,pr_palm)       = sums(nzb+14,pr_palm)       / &    ! surface temperature
2022                                    ngp_2dh(sr)
2023!--    eges, e*
2024       sums(nzb+4:nzb+5,pr_palm)  = sums(nzb+4:nzb+5,pr_palm)  / &
2025                                    ngp_3d(sr)
2026!--    Old and new divergence
2027       sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
2028                                    ngp_3d_inner(sr)
2029
2030!--    User-defined profiles
2031       IF ( max_pr_user > 0 )  THEN
2032          DO  k = nzb, nzt+1
2033             sums(k,pr_palm+1:pr_palm+max_pr_user) = &
2034                                    sums(k,pr_palm+1:pr_palm+max_pr_user) / &
2035                                    ngp_2dh_s_inner(k,sr)
2036          ENDDO
2037       ENDIF
2038
2039       IF ( air_chemistry ) THEN
2040          IF ( max_pr_cs > 0 )  THEN                 
2041             DO k = nzb, nzt+1
2042                sums(k, pr_palm+1:pr_palm+max_pr_user+max_pr_cs) = &
2043                                 sums(k, pr_palm+1:pr_palm+max_pr_user+max_pr_cs) / &
2044                                 ngp_2dh_s_inner(k,sr)
2045             ENDDO
2046          ENDIF 
2047       ENDIF
2048
2049!
2050!--    Collect horizontal average in hom.
2051!--    Compute deduced averages (e.g. total heat flux)
2052       hom(:,1,3,sr)  = sums(:,3)      ! w
2053       hom(:,1,8,sr)  = sums(:,8)      ! e     profiles 5-7 are initial profiles
2054       hom(:,1,9,sr)  = sums(:,9)      ! km
2055       hom(:,1,10,sr) = sums(:,10)     ! kh
2056       hom(:,1,11,sr) = sums(:,11)     ! l
2057       hom(:,1,12,sr) = sums(:,12)     ! w"u"
2058       hom(:,1,13,sr) = sums(:,13)     ! w*u*
2059       hom(:,1,14,sr) = sums(:,14)     ! w"v"
2060       hom(:,1,15,sr) = sums(:,15)     ! w*v*
2061       hom(:,1,16,sr) = sums(:,16)     ! w"pt"
2062       hom(:,1,17,sr) = sums(:,17)     ! w*pt*
2063       hom(:,1,18,sr) = sums(:,16) + sums(:,17)    ! wpt
2064       hom(:,1,19,sr) = sums(:,12) + sums(:,13)    ! wu
2065       hom(:,1,20,sr) = sums(:,14) + sums(:,15)    ! wv
2066       hom(:,1,21,sr) = sums(:,21)     ! w*pt*BC
2067       hom(:,1,22,sr) = sums(:,16) + sums(:,21)    ! wptBC
2068                                       ! profile 24 is initial profile (sa)
2069                                       ! profiles 25-29 left empty for initial
2070                                       ! profiles
2071       hom(:,1,30,sr) = sums(:,30)     ! u*2
2072       hom(:,1,31,sr) = sums(:,31)     ! v*2
2073       hom(:,1,32,sr) = sums(:,32)     ! w*2
2074       hom(:,1,33,sr) = sums(:,33)     ! pt*2
2075       hom(:,1,34,sr) = sums(:,34)     ! e*
2076       hom(:,1,35,sr) = sums(:,35)     ! w*2pt*
2077       hom(:,1,36,sr) = sums(:,36)     ! w*pt*2
2078       hom(:,1,37,sr) = sums(:,37)     ! w*e*
2079       hom(:,1,38,sr) = sums(:,38)     ! w*3
2080       hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20_wp )**1.5_wp   ! Sw
2081       hom(:,1,40,sr) = sums(:,40)     ! p
2082       hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
2083       hom(:,1,46,sr) = sums(:,46)     ! w*vpt*       
2084       hom(:,1,47,sr) = sums(:,45) + sums(:,46)    ! wvpt
2085       hom(:,1,48,sr) = sums(:,48)     ! w"q" (w"qv")
2086       hom(:,1,49,sr) = sums(:,49)     ! w*q* (w*qv*)
2087       hom(:,1,50,sr) = sums(:,48) + sums(:,49)    ! wq (wqv)
2088       hom(:,1,51,sr) = sums(:,51)     ! w"qv"
2089       hom(:,1,52,sr) = sums(:,52)     ! w*qv*       
2090       hom(:,1,53,sr) = sums(:,52) + sums(:,51)    ! wq (wqv)
2091       hom(:,1,54,sr) = sums(:,54)     ! ql
2092       hom(:,1,55,sr) = sums(:,55)     ! w*u*u*/dz
2093       hom(:,1,56,sr) = sums(:,56)     ! w*p*/dz
2094       hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho_ocean )/dz
2095       hom(:,1,58,sr) = sums(:,58)     ! u"pt"
2096       hom(:,1,59,sr) = sums(:,59)     ! u*pt*
2097       hom(:,1,60,sr) = sums(:,58) + sums(:,59)    ! upt_t
2098       hom(:,1,61,sr) = sums(:,61)     ! v"pt"
2099       hom(:,1,62,sr) = sums(:,62)     ! v*pt*
2100       hom(:,1,63,sr) = sums(:,61) + sums(:,62)    ! vpt_t
2101       hom(:,1,64,sr) = sums(:,64)     ! rho_ocean
2102       hom(:,1,65,sr) = sums(:,65)     ! w"sa"
2103       hom(:,1,66,sr) = sums(:,66)     ! w*sa*
2104       hom(:,1,67,sr) = sums(:,65) + sums(:,66)    ! wsa
2105       hom(:,1,68,sr) = sums(:,68)     ! w*p*
2106       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho_ocean
2107       hom(:,1,70,sr) = sums(:,70)     ! q*2
2108       hom(:,1,71,sr) = sums(:,71)     ! prho
2109       hom(:,1,72,sr) = hyp * 1E-2_wp  ! hyp in hPa
2110       hom(:,1,123,sr) = sums(:,123)   ! nc
2111       hom(:,1,73,sr) = sums(:,73)     ! nr
2112       hom(:,1,74,sr) = sums(:,74)     ! qr
2113       hom(:,1,75,sr) = sums(:,75)     ! qc
2114       hom(:,1,76,sr) = sums(:,76)     ! prr (precipitation rate)
2115                                       ! 77 is initial density profile
2116       hom(:,1,78,sr) = ug             ! ug
2117       hom(:,1,79,sr) = vg             ! vg
2118       hom(:,1,80,sr) = w_subs         ! w_subs
2119
2120       IF ( large_scale_forcing )  THEN
2121          hom(:,1,81,sr) = sums_ls_l(:,0)          ! td_lsa_lpt
2122          hom(:,1,82,sr) = sums_ls_l(:,1)          ! td_lsa_q
2123          IF ( use_subsidence_tendencies )  THEN
2124             hom(:,1,83,sr) = sums_ls_l(:,2)       ! td_sub_lpt
2125             hom(:,1,84,sr) = sums_ls_l(:,3)       ! td_sub_q
2126          ELSE
2127             hom(:,1,83,sr) = sums(:,83)           ! td_sub_lpt
2128             hom(:,1,84,sr) = sums(:,84)           ! td_sub_q
2129          ENDIF
2130          hom(:,1,85,sr) = sums(:,85)              ! td_nud_lpt
2131          hom(:,1,86,sr) = sums(:,86)              ! td_nud_q
2132          hom(:,1,87,sr) = sums(:,87)              ! td_nud_u
2133          hom(:,1,88,sr) = sums(:,88)              ! td_nud_v
2134       ENDIF
2135
2136       IF ( land_surface )  THEN
2137          hom(:,1,89,sr) = sums(:,89)              ! t_soil
2138                                                   ! 90 is initial t_soil profile
2139          hom(:,1,91,sr) = sums(:,91)              ! m_soil
2140                                                   ! 92 is initial m_soil profile
2141          hom(:,1,93,sr)  = sums(:,93)             ! ghf
2142          hom(:,1,94,sr)  = sums(:,94)             ! qsws_liq
2143          hom(:,1,95,sr)  = sums(:,95)             ! qsws_soil
2144          hom(:,1,96,sr)  = sums(:,96)             ! qsws_veg
2145          hom(:,1,97,sr)  = sums(:,97)             ! r_a
2146          hom(:,1,98,sr) = sums(:,98)              ! r_s
2147
2148       ENDIF
2149
2150       IF ( radiation )  THEN
2151          hom(:,1,99,sr) = sums(:,99)            ! rad_net
2152          hom(:,1,100,sr) = sums(:,100)            ! rad_lw_in
2153          hom(:,1,101,sr) = sums(:,101)            ! rad_lw_out
2154          hom(:,1,102,sr) = sums(:,102)            ! rad_sw_in
2155          hom(:,1,103,sr) = sums(:,103)            ! rad_sw_out
2156
2157          IF ( radiation_scheme == 'rrtmg' )  THEN
2158             hom(:,1,104,sr) = sums(:,104)            ! rad_lw_cs_hr
2159             hom(:,1,105,sr) = sums(:,105)            ! rad_lw_hr
2160             hom(:,1,106,sr) = sums(:,106)            ! rad_sw_cs_hr
2161             hom(:,1,107,sr) = sums(:,107)            ! rad_sw_hr
2162
2163             hom(:,1,108,sr) = sums(:,108)            ! rrtm_aldif
2164             hom(:,1,109,sr) = sums(:,109)            ! rrtm_aldir
2165             hom(:,1,110,sr) = sums(:,110)            ! rrtm_asdif
2166             hom(:,1,111,sr) = sums(:,111)            ! rrtm_asdir
2167          ENDIF
2168       ENDIF
2169
2170       hom(:,1,112,sr) = sums(:,112)            !: L
2171
2172       IF ( passive_scalar )  THEN
2173          hom(:,1,117,sr) = sums(:,117)     ! w"s"
2174          hom(:,1,114,sr) = sums(:,114)     ! w*s*
2175          hom(:,1,118,sr) = sums(:,117) + sums(:,114)    ! ws
2176          hom(:,1,116,sr) = sums(:,116)     ! s*2
2177       ENDIF
2178
2179       hom(:,1,119,sr) = rho_air       ! rho_air in Kg/m^3
2180       hom(:,1,120,sr) = rho_air_zw    ! rho_air_zw in Kg/m^3
2181
2182       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
2183                                       ! u*, w'u', w'v', t* (in last profile)
2184
2185       IF ( max_pr_user > 0 )  THEN    ! user-defined profiles
2186          hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
2187                               sums(:,pr_palm+1:pr_palm+max_pr_user)
2188       ENDIF
2189
2190       IF ( air_chemistry )  THEN
2191          IF ( max_pr_cs > 0 )  THEN    ! chem_spcs profiles     
2192             hom(:, 1, pr_palm+max_pr_user+1:pr_palm + max_pr_user+max_pr_cs, sr) = &
2193                               sums(:, pr_palm+max_pr_user+1:pr_palm+max_pr_user+max_pr_cs)
2194          ENDIF
2195       ENDIF
2196!
2197!--    Determine the boundary layer height using two different schemes.
2198!--    First scheme: Starting from the Earth's (Ocean's) surface, look for the
2199!--    first relative minimum (maximum) of the total heat flux.
2200!--    The corresponding height is assumed as the boundary layer height, if it
2201!--    is less than 1.5 times the height where the heat flux becomes negative
2202!--    (positive) for the first time. Attention: the resolved vertical sensible
2203!--    heat flux (hom(:,1,17,sr) = w*pt*) is not known at the beginning because
2204!--    the calculation happens in advec_s_ws which is called after
2205!--    flow_statistics. Therefore z_i is directly taken from restart data at
2206!--    the beginning of restart runs. 
2207       IF ( TRIM( initializing_actions ) /= 'read_restart_data' .OR.           &
2208            simulated_time_at_begin /= simulated_time ) THEN
2209
2210          z_i(1) = 0.0_wp
2211          first = .TRUE.
2212
2213          IF ( ocean_mode )  THEN
2214             DO  k = nzt, nzb+1, -1
2215                IF ( first  .AND.  hom(k,1,18,sr) < -1.0E-8_wp )  THEN
2216                   first = .FALSE.
2217                   height = zw(k)
2218                ENDIF
2219                IF ( hom(k,1,18,sr) < -1.0E-8_wp  .AND.                        &
2220                     hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
2221                   IF ( zw(k) < 1.5_wp * height )  THEN
2222                      z_i(1) = zw(k)
2223                   ELSE
2224                      z_i(1) = height
2225                   ENDIF
2226                   EXIT
2227                ENDIF
2228             ENDDO
2229          ELSE
2230             DO  k = nzb, nzt-1
2231                IF ( first  .AND.  hom(k,1,18,sr) < -1.0E-8_wp )  THEN
2232                   first = .FALSE.
2233                   height = zw(k)
2234                ENDIF
2235                IF ( hom(k,1,18,sr) < -1.0E-8_wp  .AND.                        &
2236                     hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
2237                   IF ( zw(k) < 1.5_wp * height )  THEN
2238                      z_i(1) = zw(k)
2239                   ELSE
2240                      z_i(1) = height
2241                   ENDIF
2242                   EXIT
2243                ENDIF
2244             ENDDO
2245          ENDIF
2246
2247!
2248!--       Second scheme: Gradient scheme from Sullivan et al. (1998), modified
2249!--       by Uhlenbrock(2006). The boundary layer height is the height with the
2250!--       maximal local temperature gradient: starting from the second (the
2251!--       last but one) vertical gridpoint, the local gradient must be at least
2252!--       0.2K/100m and greater than the next four gradients.
2253!--       WARNING: The threshold value of 0.2K/100m must be adjusted for the
2254!--       ocean case!
2255          z_i(2) = 0.0_wp
2256          DO  k = nzb+1, nzt+1
2257             dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
2258          ENDDO
2259          dptdz_threshold = 0.2_wp / 100.0_wp
2260
2261          IF ( ocean_mode )  THEN
2262             DO  k = nzt+1, nzb+5, -1
2263                IF ( dptdz(k) > dptdz_threshold  .AND.                         &
2264                     dptdz(k) > dptdz(k-1)  .AND.  dptdz(k) > dptdz(k-2)  .AND.&
2265                     dptdz(k) > dptdz(k-3)  .AND.  dptdz(k) > dptdz(k-4) )  THEN
2266                   z_i(2) = zw(k-1)
2267                   EXIT
2268                ENDIF
2269             ENDDO
2270          ELSE
2271             DO  k = nzb+1, nzt-3
2272                IF ( dptdz(k) > dptdz_threshold  .AND.                         &
2273                     dptdz(k) > dptdz(k+1)  .AND.  dptdz(k) > dptdz(k+2)  .AND.&
2274                     dptdz(k) > dptdz(k+3)  .AND.  dptdz(k) > dptdz(k+4) )  THEN
2275                   z_i(2) = zw(k-1)
2276                   EXIT
2277                ENDIF
2278             ENDDO
2279          ENDIF
2280
2281       ENDIF
2282
2283       hom(nzb+6,1,pr_palm,sr) = z_i(1)
2284       hom(nzb+7,1,pr_palm,sr) = z_i(2)
2285
2286!
2287!--    Determine vertical index which is nearest to the mean surface level
2288!--    height of the respective statistic region
2289       DO  k = nzb, nzt
2290          IF ( zw(k) >= mean_surface_level_height(sr) )  THEN
2291             k_surface_level = k
2292             EXIT
2293          ENDIF
2294       ENDDO
2295
2296!
2297!--    Computation of both the characteristic vertical velocity and
2298!--    the characteristic convective boundary layer temperature.
2299!--    The inversion height entering into the equation is defined with respect
2300!--    to the mean surface level height of the respective statistic region.
2301!--    The horizontal average at surface level index + 1 is input for the
2302!--    average temperature.
2303       IF ( hom(k_surface_level,1,18,sr) > 1.0E-8_wp  .AND.  z_i(1) /= 0.0_wp )&
2304       THEN
2305          hom(nzb+8,1,pr_palm,sr) =                                            &
2306             ( g / hom(k_surface_level+1,1,4,sr) *                             &
2307             ( hom(k_surface_level,1,18,sr) /                                  &
2308             ( heatflux_output_conversion(nzb) * rho_air(nzb) ) )              &
2309             * ABS( z_i(1) - mean_surface_level_height(sr) ) )**0.333333333_wp
2310       ELSE
2311          hom(nzb+8,1,pr_palm,sr)  = 0.0_wp
2312       ENDIF
2313
2314!
2315!--    Collect the time series quantities. Please note, timeseries quantities
2316!--    which are collected from horizontally averaged profiles, e.g. wpt
2317!--    or pt(zp), are treated specially. In case of elevated model surfaces,
2318!--    index nzb+1 might be within topography and data will be zero. Therefore,
2319!--    take value for the first atmosphere index, which is topo_min_level+1.
2320       ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)        ! E
2321       ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)        ! E*
2322       ts_value(3,sr) = dt_3d
2323       ts_value(4,sr) = hom(nzb,1,pr_palm,sr)          ! u*
2324       ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr)        ! th*
2325       ts_value(6,sr) = u_max
2326       ts_value(7,sr) = v_max
2327       ts_value(8,sr) = w_max
2328       ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr)       ! new divergence
2329       ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr)       ! old Divergence
2330       ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr)       ! z_i(1)
2331       ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr)       ! z_i(2)
2332       ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr)       ! w*
2333       ts_value(14,sr) = hom(nzb,1,16,sr)              ! w'pt'   at k=0
2334       ts_value(15,sr) = hom(topo_min_level+1,1,16,sr) ! w'pt'   at k=1
2335       ts_value(16,sr) = hom(topo_min_level+1,1,18,sr) ! wpt     at k=1
2336       ts_value(17,sr) = hom(nzb+14,1,pr_palm,sr)      ! pt(0)
2337       ts_value(18,sr) = hom(topo_min_level+1,1,4,sr)  ! pt(zp)
2338       ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr)       ! u'w'    at k=0
2339       ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)       ! v'w'    at k=0
2340       ts_value(21,sr) = hom(nzb,1,48,sr)              ! w"q"    at k=0
2341
2342       IF ( .NOT. neutral )  THEN
2343          ts_value(22,sr) = hom(nzb,1,112,sr)          ! L
2344       ELSE
2345          ts_value(22,sr) = 1.0E10_wp
2346       ENDIF
2347
2348       ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
2349
2350       IF ( passive_scalar )  THEN
2351          ts_value(24,sr) = hom(nzb+13,1,117,sr)       ! w"s" ( to do ! )
2352          ts_value(25,sr) = hom(nzb+13,1,pr_palm,sr)   ! s*
2353       ENDIF
2354
2355!
2356!--    Collect land surface model timeseries
2357       IF ( land_surface )  THEN
2358          ts_value(dots_soil  ,sr) = hom(nzb,1,93,sr)           ! ghf
2359          ts_value(dots_soil+1,sr) = hom(nzb,1,94,sr)           ! qsws_liq
2360          ts_value(dots_soil+2,sr) = hom(nzb,1,95,sr)           ! qsws_soil
2361          ts_value(dots_soil+3,sr) = hom(nzb,1,96,sr)           ! qsws_veg
2362          ts_value(dots_soil+4,sr) = hom(nzb,1,97,sr)           ! r_a
2363          ts_value(dots_soil+5,sr) = hom(nzb,1,98,sr)           ! r_s
2364       ENDIF
2365!
2366!--    Collect radiation model timeseries
2367       IF ( radiation )  THEN
2368          ts_value(dots_rad,sr)   = hom(nzb,1,99,sr)           ! rad_net
2369          ts_value(dots_rad+1,sr) = hom(nzb,1,100,sr)          ! rad_lw_in
2370          ts_value(dots_rad+2,sr) = hom(nzb,1,101,sr)          ! rad_lw_out
2371          ts_value(dots_rad+3,sr) = hom(nzb,1,102,sr)          ! rad_sw_in
2372          ts_value(dots_rad+4,sr) = hom(nzb,1,103,sr)          ! rad_sw_out
2373
2374          IF ( radiation_scheme == 'rrtmg' )  THEN
2375             ts_value(dots_rad+5,sr) = hom(nzb,1,108,sr)          ! rrtm_aldif
2376             ts_value(dots_rad+6,sr) = hom(nzb,1,109,sr)          ! rrtm_aldir
2377             ts_value(dots_rad+7,sr) = hom(nzb,1,110,sr)          ! rrtm_asdif
2378             ts_value(dots_rad+8,sr) = hom(nzb,1,111,sr)          ! rrtm_asdir
2379          ENDIF
2380
2381       ENDIF
2382
2383!
2384!--    Calculate additional statistics provided by other modules
2385       CALL module_interface_statistics( 'time_series', sr, 0, dots_max )
2386
2387    ENDDO    ! loop of the subregions
2388
2389!
2390!-- If required, sum up horizontal averages for subsequent time averaging.
2391!-- Do not sum, if flow statistics is called before the first initial time step.
2392    IF ( do_sum  .AND.  simulated_time /= 0.0_wp )  THEN
2393       IF ( average_count_pr == 0 )  hom_sum = 0.0_wp
2394       hom_sum = hom_sum + hom(:,1,:,:)
2395       average_count_pr = average_count_pr + 1
2396       do_sum = .FALSE.
2397    ENDIF
2398
2399!
2400!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
2401!-- This flag is reset after each time step in time_integration.
2402    flow_statistics_called = .TRUE.
2403
2404    CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
2405
2406
2407 END SUBROUTINE flow_statistics
Note: See TracBrowser for help on using the repository browser.