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

Last change on this file since 3828 was 3828, checked in by raasch, 5 years ago

further gfortran warnings activated on testserver, unused variables removed

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