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

Last change on this file since 1901 was 1854, checked in by maronga, 8 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 147.0 KB
Line 
1!> @file flow_statistics.f90
2!--------------------------------------------------------------------------------!
3! This file is part of PALM.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms
6! of the GNU General Public License as published by the Free Software Foundation,
7! either version 3 of the License, or (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with
14! PALM. If not, see <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2016 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------!
18!
19! Current revisions:
20! -----------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: flow_statistics.f90 1854 2016-04-11 09:03:15Z raasch $
26!
27! 1853 2016-04-11 09:00:35Z maronga
28! Adjusted for use with radiation_scheme = constant
29!
30! 1849 2016-04-08 11:33:18Z hoffmann
31! prr moved to arrays_3d
32!
33! 1822 2016-04-07 07:49:42Z hoffmann
34! Output of bulk microphysics simplified.
35!
36! 1815 2016-04-06 13:49:59Z raasch
37! cpp-directives for intel openmp bug removed
38!
39! 1783 2016-03-06 18:36:17Z raasch
40! +module netcdf_interface
41!
42! 1747 2016-02-08 12:25:53Z raasch
43! small bugfixes for accelerator version
44!
45! 1738 2015-12-18 13:56:05Z raasch
46! bugfixes for calculations in statistical regions which do not contain grid
47! points in the lowest vertical levels, mean surface level height considered
48! in the calculation of the characteristic vertical velocity,
49! old upstream parts removed
50!
51! 1709 2015-11-04 14:47:01Z maronga
52! Updated output of Obukhov length
53!
54! 1691 2015-10-26 16:17:44Z maronga
55! Revised calculation of Obukhov length. Added output of radiative heating >
56! rates for RRTMG.
57!
58! 1682 2015-10-07 23:56:08Z knoop
59! Code annotations made doxygen readable
60!
61! 1658 2015-09-18 10:52:53Z raasch
62! bugfix: temporary reduction variables in the openacc branch are now
63! initialized to zero
64!
65! 1654 2015-09-17 09:20:17Z raasch
66! FORTRAN bugfix of r1652
67!
68! 1652 2015-09-17 08:12:24Z raasch
69! bugfix in calculation of energy production by turbulent transport of TKE
70!
71! 1593 2015-05-16 13:58:02Z raasch
72! FORTRAN errors removed from openacc branch
73!
74! 1585 2015-04-30 07:05:52Z maronga
75! Added output of timeseries and profiles for RRTMG
76!
77! 1571 2015-03-12 16:12:49Z maronga
78! Bugfix: output of rad_net and rad_sw_in
79!
80! 1567 2015-03-10 17:57:55Z suehring
81! Reverse modifications made for monotonic limiter.
82!
83! 1557 2015-03-05 16:43:04Z suehring
84! Adjustments for monotonic limiter
85!
86! 1555 2015-03-04 17:44:27Z maronga
87! Added output of r_a and r_s.
88!
89! 1551 2015-03-03 14:18:16Z maronga
90! Added suppport for land surface model and radiation model output.
91!
92! 1498 2014-12-03 14:09:51Z suehring
93! Comments added
94!
95! 1482 2014-10-18 12:34:45Z raasch
96! missing ngp_sums_ls added in accelerator version
97!
98! 1450 2014-08-21 07:31:51Z heinze
99! bugfix: calculate fac only for simulated_time >= 0.0
100!
101! 1396 2014-05-06 13:37:41Z raasch
102! bugfix: "copyin" replaced by "update device" in openacc-branch
103!
104! 1386 2014-05-05 13:55:30Z boeske
105! bugfix: simulated time before the last timestep is needed for the correct
106! calculation of the profiles of large scale forcing tendencies
107!
108! 1382 2014-04-30 12:15:41Z boeske
109! Renamed variables which store large scale forcing tendencies
110! pt_lsa -> td_lsa_lpt, pt_subs -> td_sub_lpt,
111! q_lsa  -> td_lsa_q,   q_subs  -> td_sub_q,
112! added Neumann boundary conditions for profile data output of large scale
113! advection and subsidence terms at nzt+1
114!
115! 1374 2014-04-25 12:55:07Z raasch
116! bugfix: syntax errors removed from openacc-branch
117! missing variables added to ONLY-lists
118!
119! 1365 2014-04-22 15:03:56Z boeske
120! Output of large scale advection, large scale subsidence and nudging tendencies
121! +sums_ls_l, ngp_sums_ls, use_subsidence_tendencies
122!
123! 1353 2014-04-08 15:21:23Z heinze
124! REAL constants provided with KIND-attribute
125!
126! 1322 2014-03-20 16:38:49Z raasch
127! REAL constants defined as wp-kind
128!
129! 1320 2014-03-20 08:40:49Z raasch
130! ONLY-attribute added to USE-statements,
131! kind-parameters added to all INTEGER and REAL declaration statements,
132! kinds are defined in new module kinds,
133! revision history before 2012 removed,
134! comment fields (!:) to be used for variable explanations added to
135! all variable declaration statements
136!
137! 1299 2014-03-06 13:15:21Z heinze
138! Output of large scale vertical velocity w_subs
139!
140! 1257 2013-11-08 15:18:40Z raasch
141! openacc "end parallel" replaced by "end parallel loop"
142!
143! 1241 2013-10-30 11:36:58Z heinze
144! Output of ug and vg
145!
146! 1221 2013-09-10 08:59:13Z raasch
147! ported for openACC in separate #else branch
148!
149! 1179 2013-06-14 05:57:58Z raasch
150! comment for profile 77 added
151!
152! 1115 2013-03-26 18:16:16Z hoffmann
153! ql is calculated by calc_liquid_water_content
154!
155! 1111 2013-03-08 23:54:10Z raasch
156! openACC directive added
157!
158! 1053 2012-11-13 17:11:03Z hoffmann
159! additions for two-moment cloud physics scheme:
160! +nr, qr, qc, prr
161!
162! 1036 2012-10-22 13:43:42Z raasch
163! code put under GPL (PALM 3.9)
164!
165! 1007 2012-09-19 14:30:36Z franke
166! Calculation of buoyancy flux for humidity in case of WS-scheme is now using
167! turbulent fluxes of WS-scheme
168! Bugfix: Calculation of subgridscale buoyancy flux for humidity and cloud
169! droplets at nzb and nzt added
170!
171! 801 2012-01-10 17:30:36Z suehring
172! Calculation of turbulent fluxes in advec_ws is now thread-safe.
173!
174! Revision 1.1  1997/08/11 06:15:17  raasch
175! Initial revision
176!
177!
178! Description:
179! ------------
180!> Compute average profiles and further average flow quantities for the different
181!> user-defined (sub-)regions. The region indexed 0 is the total model domain.
182!>
183!> @note For simplicity, nzb_s_inner and nzb_diff_s_inner are being used as a
184!>       lower vertical index for k-loops for all variables, although strictly
185!>       speaking the k-loops would have to be split up according to the staggered
186!>       grid. However, this implies no error since staggered velocity components
187!>       are zero at the walls and inside buildings.
188!------------------------------------------------------------------------------!
189#if ! defined( __openacc )
190 SUBROUTINE flow_statistics
191 
192
193    USE arrays_3d,                                                             &
194        ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, ol, p, prho, prr, pt, q, qc, ql,&
195               qr, qs, qsws, qswst, rho, sa, saswsb, saswst, shf, td_lsa_lpt,  &
196               td_lsa_q, td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us,& 
197               usws, uswst, vsws, v, vg, vpt, vswst, w, w_subs, zw
198       
199    USE cloud_parameters,                                                      &
200        ONLY:   l_d_cp, pt_d_t
201       
202    USE control_parameters,                                                    &
203        ONLY:   average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
204                dt_3d, g, humidity, kappa, large_scale_forcing,                &
205                large_scale_subsidence, max_pr_user, message_string, neutral,  &
206                microphysics_seifert, ocean, passive_scalar, simulated_time,   &
207                use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
208                ws_scheme_mom, ws_scheme_sca
209       
210    USE cpulog,                                                                &
211        ONLY:   cpu_log, log_point
212       
213    USE grid_variables,                                                        &
214        ONLY:   ddx, ddy
215       
216    USE indices,                                                               &
217        ONLY:   ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums,      &
218                ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzb_diff_s_inner,        &
219                nzb_s_inner, nzt, nzt_diff
220       
221    USE kinds
222   
223    USE land_surface_model_mod,                                                &
224        ONLY:   ghf_eb, land_surface, m_soil, nzb_soil, nzt_soil,              &
225                qsws_eb, qsws_liq_eb, qsws_soil_eb, qsws_veg_eb, r_a, r_s,     &
226                shf_eb, t_soil
227
228    USE netcdf_interface,                                                      &
229        ONLY:  dots_rad, dots_soil
230
231    USE pegrid
232
233    USE radiation_model_mod,                                                   &
234        ONLY:  radiation, radiation_scheme, rad_net,                           &
235               rad_lw_in, rad_lw_out, rad_lw_cs_hr, rad_lw_hr,                 &
236               rad_sw_in, rad_sw_out, rad_sw_cs_hr, rad_sw_hr
237
238#if defined ( __rrtmg )
239    USE radiation_model_mod,                                                   &
240        ONLY:  rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
241#endif
242 
243    USE statistics
244
245
246    IMPLICIT NONE
247
248    INTEGER(iwp) ::  i                   !<
249    INTEGER(iwp) ::  j                   !<
250    INTEGER(iwp) ::  k                   !<
251    INTEGER(iwp) ::  k_surface_level     !<
252    INTEGER(iwp) ::  nt                  !<
253    INTEGER(iwp) ::  omp_get_thread_num  !<
254    INTEGER(iwp) ::  sr                  !<
255    INTEGER(iwp) ::  tn                  !<
256   
257    LOGICAL ::  first  !<
258   
259    REAL(wp) ::  dptdz_threshold  !<
260    REAL(wp) ::  fac              !<
261    REAL(wp) ::  height           !<
262    REAL(wp) ::  pts              !<
263    REAL(wp) ::  sums_l_eper      !<
264    REAL(wp) ::  sums_l_etot      !<
265    REAL(wp) ::  ust              !<
266    REAL(wp) ::  ust2             !<
267    REAL(wp) ::  u2               !<
268    REAL(wp) ::  vst              !<
269    REAL(wp) ::  vst2             !<
270    REAL(wp) ::  v2               !<
271    REAL(wp) ::  w2               !<
272    REAL(wp) ::  z_i(2)           !<
273   
274    REAL(wp) ::  dptdz(nzb+1:nzt+1)    !<
275    REAL(wp) ::  sums_ll(nzb:nzt+1,2)  !<
276
277    CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
278
279    !$acc update host( km, kh, e, ol, pt, qs, qsws, shf, ts, u, usws, v, vsws, w )
280
281!
282!-- To be on the safe side, check whether flow_statistics has already been
283!-- called once after the current time step
284    IF ( flow_statistics_called )  THEN
285
286       message_string = 'flow_statistics is called two times within one ' // &
287                        'timestep'
288       CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
289
290    ENDIF
291
292!
293!-- Compute statistics for each (sub-)region
294    DO  sr = 0, statistic_regions
295
296!
297!--    Initialize (local) summation array
298       sums_l = 0.0_wp
299
300!
301!--    Store sums that have been computed in other subroutines in summation
302!--    array
303       sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
304!--    WARNING: next line still has to be adjusted for OpenMP
305       sums_l(:,21,0) = sums_wsts_bc_l(:,sr)  ! heat flux from advec_s_bc
306       sums_l(nzb+9,pr_palm,0)  = sums_divold_l(sr)  ! old divergence from pres
307       sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
308
309!
310!--    When calcuating horizontally-averaged total (resolved- plus subgrid-
311!--    scale) vertical fluxes and velocity variances by using commonly-
312!--    applied Reynolds-based methods ( e.g. <w'pt'> = (w-<w>)*(pt-<pt>) )
313!--    in combination with the 5th order advection scheme, pronounced
314!--    artificial kinks could be observed in the vertical profiles near the
315!--    surface. Please note: these kinks were not related to the model truth,
316!--    i.e. these kinks are just related to an evaluation problem.   
317!--    In order avoid these kinks, vertical fluxes and horizontal as well
318!--    vertical velocity variances are calculated directly within the advection
319!--    routines, according to the numerical discretization, to evaluate the
320!--    statistical quantities as they will appear within the prognostic
321!--    equations.
322!--    Copy the turbulent quantities, evaluated in the advection routines to
323!--    the local array sums_l() for further computations.
324       IF ( ws_scheme_mom .AND. sr == 0 )  THEN
325
326!
327!--       According to the Neumann bc for the horizontal velocity components,
328!--       the corresponding fluxes has to satisfiy the same bc.
329          IF ( ocean )  THEN
330             sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
331             sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)
332          ENDIF
333
334          DO  i = 0, threads_per_task-1
335!
336!--          Swap the turbulent quantities evaluated in advec_ws.
337             sums_l(:,13,i) = sums_wsus_ws_l(:,i)       ! w*u*
338             sums_l(:,15,i) = sums_wsvs_ws_l(:,i)       ! w*v*
339             sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2
340             sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2
341             sums_l(:,32,i) = sums_ws2_ws_l(:,i)        ! w*2
342             sums_l(:,34,i) = sums_l(:,34,i) + 0.5_wp *                        & 
343                              ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +      &
344                                sums_ws2_ws_l(:,i) )    ! e*
345             DO  k = nzb, nzt
346                sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5_wp * ( &
347                                                      sums_us2_ws_l(k,i) +     &
348                                                      sums_vs2_ws_l(k,i) +     &
349                                                      sums_ws2_ws_l(k,i) )
350             ENDDO
351          ENDDO
352
353       ENDIF
354
355       IF ( ws_scheme_sca .AND. sr == 0 )  THEN
356
357          DO  i = 0, threads_per_task-1
358             sums_l(:,17,i) = sums_wspts_ws_l(:,i)      ! w*pt* from advec_s_ws
359             IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
360             IF ( humidity .OR. passive_scalar ) sums_l(:,49,i) =              &
361                                                   sums_wsqs_ws_l(:,i) !w*q*
362          ENDDO
363
364       ENDIF
365!
366!--    Horizontally averaged profiles of horizontal velocities and temperature.
367!--    They must have been computed before, because they are already required
368!--    for other horizontal averages.
369       tn = 0
370
371       !$OMP PARALLEL PRIVATE( i, j, k, tn )
372!$     tn = omp_get_thread_num()
373
374       !$OMP DO
375       DO  i = nxl, nxr
376          DO  j =  nys, nyn
377             DO  k = nzb_s_inner(j,i), nzt+1
378                sums_l(k,1,tn)  = sums_l(k,1,tn)  + u(k,j,i)  * rmask(j,i,sr)
379                sums_l(k,2,tn)  = sums_l(k,2,tn)  + v(k,j,i)  * rmask(j,i,sr)
380                sums_l(k,4,tn)  = sums_l(k,4,tn)  + pt(k,j,i) * rmask(j,i,sr)
381             ENDDO
382          ENDDO
383       ENDDO
384
385!
386!--    Horizontally averaged profile of salinity
387       IF ( ocean )  THEN
388          !$OMP DO
389          DO  i = nxl, nxr
390             DO  j =  nys, nyn
391                DO  k = nzb_s_inner(j,i), nzt+1
392                   sums_l(k,23,tn)  = sums_l(k,23,tn) + &
393                                      sa(k,j,i) * rmask(j,i,sr)
394                ENDDO
395             ENDDO
396          ENDDO
397       ENDIF
398
399!
400!--    Horizontally averaged profiles of virtual potential temperature,
401!--    total water content, specific humidity and liquid water potential
402!--    temperature
403       IF ( humidity )  THEN
404          !$OMP DO
405          DO  i = nxl, nxr
406             DO  j =  nys, nyn
407                DO  k = nzb_s_inner(j,i), nzt+1
408                   sums_l(k,44,tn)  = sums_l(k,44,tn) + &
409                                      vpt(k,j,i) * rmask(j,i,sr)
410                   sums_l(k,41,tn)  = sums_l(k,41,tn) + &
411                                      q(k,j,i) * rmask(j,i,sr)
412                ENDDO
413             ENDDO
414          ENDDO
415          IF ( cloud_physics )  THEN
416             !$OMP DO
417             DO  i = nxl, nxr
418                DO  j =  nys, nyn
419                   DO  k = nzb_s_inner(j,i), nzt+1
420                      sums_l(k,42,tn) = sums_l(k,42,tn) + &
421                                      ( q(k,j,i) - ql(k,j,i) ) * rmask(j,i,sr)
422                      sums_l(k,43,tn) = sums_l(k,43,tn) + ( &
423                                      pt(k,j,i) + l_d_cp*pt_d_t(k) * ql(k,j,i) &
424                                                          ) * rmask(j,i,sr)
425                   ENDDO
426                ENDDO
427             ENDDO
428          ENDIF
429       ENDIF
430
431!
432!--    Horizontally averaged profiles of passive scalar
433       IF ( passive_scalar )  THEN
434          !$OMP DO
435          DO  i = nxl, nxr
436             DO  j =  nys, nyn
437                DO  k = nzb_s_inner(j,i), nzt+1
438                   sums_l(k,41,tn)  = sums_l(k,41,tn) + q(k,j,i) * rmask(j,i,sr)
439                ENDDO
440             ENDDO
441          ENDDO
442       ENDIF
443       !$OMP END PARALLEL
444!
445!--    Summation of thread sums
446       IF ( threads_per_task > 1 )  THEN
447          DO  i = 1, threads_per_task-1
448             sums_l(:,1,0) = sums_l(:,1,0) + sums_l(:,1,i)
449             sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i)
450             sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i)
451             IF ( ocean )  THEN
452                sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)
453             ENDIF
454             IF ( humidity )  THEN
455                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
456                sums_l(:,44,0) = sums_l(:,44,0) + sums_l(:,44,i)
457                IF ( cloud_physics )  THEN
458                   sums_l(:,42,0) = sums_l(:,42,0) + sums_l(:,42,i)
459                   sums_l(:,43,0) = sums_l(:,43,0) + sums_l(:,43,i)
460                ENDIF
461             ENDIF
462             IF ( passive_scalar )  THEN
463                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
464             ENDIF
465          ENDDO
466       ENDIF
467
468#if defined( __parallel )
469!
470!--    Compute total sum from local sums
471       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
472       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL,  &
473                           MPI_SUM, comm2d, ierr )
474       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
475       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL,  &
476                           MPI_SUM, comm2d, ierr )
477       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
478       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL,  &
479                           MPI_SUM, comm2d, ierr )
480       IF ( ocean )  THEN
481          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
482          CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb,       &
483                              MPI_REAL, MPI_SUM, comm2d, ierr )
484       ENDIF
485       IF ( humidity ) THEN
486          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
487          CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb,       &
488                              MPI_REAL, MPI_SUM, comm2d, ierr )
489          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
490          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb,       &
491                              MPI_REAL, MPI_SUM, comm2d, ierr )
492          IF ( cloud_physics ) THEN
493             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
494             CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb,    &
495                                 MPI_REAL, MPI_SUM, comm2d, ierr )
496             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
497             CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb,    &
498                                 MPI_REAL, MPI_SUM, comm2d, ierr )
499          ENDIF
500       ENDIF
501
502       IF ( passive_scalar )  THEN
503          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
504          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb,       &
505                              MPI_REAL, MPI_SUM, comm2d, ierr )
506       ENDIF
507#else
508       sums(:,1) = sums_l(:,1,0)
509       sums(:,2) = sums_l(:,2,0)
510       sums(:,4) = sums_l(:,4,0)
511       IF ( ocean )  sums(:,23) = sums_l(:,23,0)
512       IF ( humidity ) THEN
513          sums(:,44) = sums_l(:,44,0)
514          sums(:,41) = sums_l(:,41,0)
515          IF ( cloud_physics ) THEN
516             sums(:,42) = sums_l(:,42,0)
517             sums(:,43) = sums_l(:,43,0)
518          ENDIF
519       ENDIF
520       IF ( passive_scalar )  sums(:,41) = sums_l(:,41,0)
521#endif
522
523!
524!--    Final values are obtained by division by the total number of grid points
525!--    used for summation. After that store profiles.
526       sums(:,1) = sums(:,1) / ngp_2dh(sr)
527       sums(:,2) = sums(:,2) / ngp_2dh(sr)
528       sums(:,4) = sums(:,4) / ngp_2dh_s_inner(:,sr)
529       hom(:,1,1,sr) = sums(:,1)             ! u
530       hom(:,1,2,sr) = sums(:,2)             ! v
531       hom(:,1,4,sr) = sums(:,4)             ! pt
532
533
534!
535!--    Salinity
536       IF ( ocean )  THEN
537          sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr)
538          hom(:,1,23,sr) = sums(:,23)             ! sa
539       ENDIF
540
541!
542!--    Humidity and cloud parameters
543       IF ( humidity ) THEN
544          sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr)
545          sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
546          hom(:,1,44,sr) = sums(:,44)             ! vpt
547          hom(:,1,41,sr) = sums(:,41)             ! qv (q)
548          IF ( cloud_physics ) THEN
549             sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr)
550             sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr)
551             hom(:,1,42,sr) = sums(:,42)             ! qv
552             hom(:,1,43,sr) = sums(:,43)             ! pt
553          ENDIF
554       ENDIF
555
556!
557!--    Passive scalar
558       IF ( passive_scalar )  hom(:,1,41,sr) = sums(:,41) /                    &
559            ngp_2dh_s_inner(:,sr)                    ! s (q)
560
561!
562!--    Horizontally averaged profiles of the remaining prognostic variables,
563!--    variances, the total and the perturbation energy (single values in last
564!--    column of sums_l) and some diagnostic quantities.
565!--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
566!--    ----  speaking the following k-loop would have to be split up and
567!--          rearranged according to the staggered grid.
568!--          However, this implies no error since staggered velocity components
569!--          are zero at the walls and inside buildings.
570       tn = 0
571       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper,             &
572       !$OMP                   sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2,  &
573       !$OMP                   w2 )
574!$     tn = omp_get_thread_num()
575
576       !$OMP DO
577       DO  i = nxl, nxr
578          DO  j =  nys, nyn
579             sums_l_etot = 0.0_wp
580             DO  k = nzb_s_inner(j,i), nzt+1
581!
582!--             Prognostic and diagnostic variables
583                sums_l(k,3,tn)  = sums_l(k,3,tn)  + w(k,j,i)  * rmask(j,i,sr)
584                sums_l(k,8,tn)  = sums_l(k,8,tn)  + e(k,j,i)  * rmask(j,i,sr)
585                sums_l(k,9,tn)  = sums_l(k,9,tn)  + km(k,j,i) * rmask(j,i,sr)
586                sums_l(k,10,tn) = sums_l(k,10,tn) + kh(k,j,i) * rmask(j,i,sr)
587                sums_l(k,40,tn) = sums_l(k,40,tn) + p(k,j,i)
588
589                sums_l(k,33,tn) = sums_l(k,33,tn) + &
590                                  ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr)
591
592                IF ( humidity )  THEN
593                   sums_l(k,70,tn) = sums_l(k,70,tn) + &
594                                  ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr)
595                ENDIF
596
597!
598!--             Higher moments
599!--             (Computation of the skewness of w further below)
600                sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i)**3 * rmask(j,i,sr)
601
602                sums_l_etot  = sums_l_etot + &
603                                        0.5_wp * ( u(k,j,i)**2 + v(k,j,i)**2 + &
604                                        w(k,j,i)**2 ) * rmask(j,i,sr)
605             ENDDO
606!
607!--          Total and perturbation energy for the total domain (being
608!--          collected in the last column of sums_l). Summation of these
609!--          quantities is seperated from the previous loop in order to
610!--          allow vectorization of that loop.
611             sums_l(nzb+4,pr_palm,tn) = sums_l(nzb+4,pr_palm,tn) + sums_l_etot
612!
613!--          2D-arrays (being collected in the last column of sums_l)
614             sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +               &
615                                        us(j,i)   * rmask(j,i,sr)
616             sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +             &
617                                        usws(j,i) * rmask(j,i,sr)
618             sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +             &
619                                        vsws(j,i) * rmask(j,i,sr)
620             sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +             &
621                                        ts(j,i)   * rmask(j,i,sr)
622             IF ( humidity )  THEN
623                sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +        &
624                                            qs(j,i)   * rmask(j,i,sr)
625             ENDIF
626          ENDDO
627       ENDDO
628
629!
630!--    Computation of statistics when ws-scheme is not used. Else these
631!--    quantities are evaluated in the advection routines.
632       IF ( .NOT. ws_scheme_mom .OR. sr /= 0 )  THEN
633          !$OMP DO
634          DO  i = nxl, nxr
635             DO  j =  nys, nyn
636                sums_l_eper = 0.0_wp
637                DO  k = nzb_s_inner(j,i), nzt+1
638                   u2   = u(k,j,i)**2
639                   v2   = v(k,j,i)**2
640                   w2   = w(k,j,i)**2
641                   ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
642                   vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
643
644                   sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr)
645                   sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr)
646                   sums_l(k,32,tn) = sums_l(k,32,tn) + w2   * rmask(j,i,sr)
647!
648!--             Perturbation energy
649
650                   sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5_wp *                &
651                                  ( ust2 + vst2 + w2 ) * rmask(j,i,sr)
652                   sums_l_eper     = sums_l_eper +                             &
653                                     0.5_wp * ( ust2+vst2+w2 ) * rmask(j,i,sr)
654
655                ENDDO
656                sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn)            &
657                     + sums_l_eper
658             ENDDO
659          ENDDO
660       ENDIF
661
662!
663!--    Horizontally averaged profiles of the vertical fluxes
664
665       !$OMP DO
666       DO  i = nxl, nxr
667          DO  j = nys, nyn
668!
669!--          Subgridscale fluxes (without Prandtl layer from k=nzb,
670!--          oterwise from k=nzb+1)
671!--          NOTE: for simplicity, nzb_diff_s_inner is used below, although
672!--          ----  strictly speaking the following k-loop would have to be
673!--                split up according to the staggered grid.
674!--                However, this implies no error since staggered velocity
675!--                components are zero at the walls and inside buildings.
676
677             DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
678!
679!--             Momentum flux w"u"
680                sums_l(k,12,tn) = sums_l(k,12,tn) - 0.25_wp * (                &
681                               km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
682                                                           ) * (               &
683                                   ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
684                                 + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
685                                                               ) * rmask(j,i,sr)
686!
687!--             Momentum flux w"v"
688                sums_l(k,14,tn) = sums_l(k,14,tn) - 0.25_wp * (                &
689                               km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
690                                                           ) * (               &
691                                   ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
692                                 + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
693                                                               ) * rmask(j,i,sr)
694!
695!--             Heat flux w"pt"
696                sums_l(k,16,tn) = sums_l(k,16,tn)                              &
697                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
698                                               * ( pt(k+1,j,i) - pt(k,j,i) )   &
699                                               * ddzu(k+1) * rmask(j,i,sr)
700
701
702!
703!--             Salinity flux w"sa"
704                IF ( ocean )  THEN
705                   sums_l(k,65,tn) = sums_l(k,65,tn)                           &
706                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
707                                               * ( sa(k+1,j,i) - sa(k,j,i) )   &
708                                               * ddzu(k+1) * rmask(j,i,sr)
709                ENDIF
710
711!
712!--             Buoyancy flux, water flux (humidity flux) w"q"
713                IF ( humidity ) THEN
714                   sums_l(k,45,tn) = sums_l(k,45,tn)                           &
715                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
716                                               * ( vpt(k+1,j,i) - vpt(k,j,i) ) &
717                                               * ddzu(k+1) * rmask(j,i,sr)
718                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
719                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
720                                               * ( q(k+1,j,i) - q(k,j,i) )     &
721                                               * ddzu(k+1) * rmask(j,i,sr)
722
723                   IF ( cloud_physics ) THEN
724                      sums_l(k,51,tn) = sums_l(k,51,tn)                        &
725                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
726                                               * ( ( q(k+1,j,i) - ql(k+1,j,i) )&
727                                                - ( q(k,j,i) - ql(k,j,i) ) )   &
728                                               * ddzu(k+1) * rmask(j,i,sr) 
729                   ENDIF
730                ENDIF
731
732!
733!--             Passive scalar flux
734                IF ( passive_scalar )  THEN
735                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
736                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
737                                               * ( q(k+1,j,i) - q(k,j,i) )     &
738                                               * ddzu(k+1) * rmask(j,i,sr)
739                ENDIF
740
741             ENDDO
742
743!
744!--          Subgridscale fluxes in the Prandtl layer
745             IF ( use_surface_fluxes )  THEN
746                sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
747                                    usws(j,i) * rmask(j,i,sr)     ! w"u"
748                sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
749                                    vsws(j,i) * rmask(j,i,sr)     ! w"v"
750                sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
751                                    shf(j,i)  * rmask(j,i,sr)     ! w"pt"
752                sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
753                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
754                sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
755                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
756                IF ( ocean )  THEN
757                   sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
758                                       saswsb(j,i) * rmask(j,i,sr)  ! w"sa"
759                ENDIF
760                IF ( humidity )  THEN
761                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
762                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
763                   sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                   &
764                                       ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) *     &
765                                       shf(j,i) + 0.61_wp * pt(nzb,j,i) *      &
766                                                  qsws(j,i) )
767                   IF ( cloud_droplets )  THEN
768                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                &
769                                         ( 1.0_wp + 0.61_wp * q(nzb,j,i) -     &
770                                           ql(nzb,j,i) ) * shf(j,i) +          &
771                                           0.61_wp * pt(nzb,j,i) * qsws(j,i) )
772                   ENDIF
773                   IF ( cloud_physics )  THEN
774!
775!--                   Formula does not work if ql(nzb) /= 0.0
776                      sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                  & 
777                                          qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
778                   ENDIF
779                ENDIF
780                IF ( passive_scalar )  THEN
781                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
782                                       qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
783                ENDIF
784             ENDIF
785
786             IF ( .NOT. neutral )  THEN
787                sums_l(nzb,114,tn) = sums_l(nzb,114,tn) +                      &
788                                    ol(j,i)  * rmask(j,i,sr) ! L
789             ENDIF
790
791
792             IF ( land_surface )  THEN
793                sums_l(nzb,93,tn)  = sums_l(nzb,93,tn) + ghf_eb(j,i)
794                sums_l(nzb,94,tn)  = sums_l(nzb,94,tn) + shf_eb(j,i)
795                sums_l(nzb,95,tn)  = sums_l(nzb,95,tn) + qsws_eb(j,i)
796                sums_l(nzb,96,tn)  = sums_l(nzb,96,tn) + qsws_liq_eb(j,i)
797                sums_l(nzb,97,tn)  = sums_l(nzb,97,tn) + qsws_soil_eb(j,i)
798                sums_l(nzb,98,tn)  = sums_l(nzb,98,tn) + qsws_veg_eb(j,i)
799                sums_l(nzb,99,tn)  = sums_l(nzb,99,tn) + r_a(j,i)
800                sums_l(nzb,100,tn) = sums_l(nzb,100,tn)+ r_s(j,i)
801             ENDIF
802
803             IF ( radiation  .AND.  radiation_scheme /= 'constant' )  THEN
804                sums_l(nzb,101,tn)  = sums_l(nzb,101,tn)  + rad_net(j,i)
805                sums_l(nzb,102,tn)  = sums_l(nzb,102,tn)  + rad_lw_in(nzb,j,i)
806                sums_l(nzb,103,tn)  = sums_l(nzb,103,tn)  + rad_lw_out(nzb,j,i)
807                sums_l(nzb,104,tn)  = sums_l(nzb,104,tn)  + rad_sw_in(nzb,j,i)
808                sums_l(nzb,105,tn)  = sums_l(nzb,105,tn)  + rad_sw_out(nzb,j,i)
809
810#if defined ( __rrtmg )
811                IF ( radiation_scheme == 'rrtmg' )  THEN
812                   sums_l(nzb,110,tn)  = sums_l(nzb,110,tn)  + rrtm_aldif(0,j,i)
813                   sums_l(nzb,111,tn)  = sums_l(nzb,111,tn)  + rrtm_aldir(0,j,i)
814                   sums_l(nzb,112,tn)  = sums_l(nzb,112,tn)  + rrtm_asdif(0,j,i)
815                   sums_l(nzb,113,tn)  = sums_l(nzb,113,tn)  + rrtm_asdir(0,j,i)
816                ENDIF
817#endif
818
819             ENDIF
820
821!
822!--          Subgridscale fluxes at the top surface
823             IF ( use_top_fluxes )  THEN
824                sums_l(nzt:nzt+1,12,tn) = sums_l(nzt:nzt+1,12,tn) + &
825                                    uswst(j,i) * rmask(j,i,sr)    ! w"u"
826                sums_l(nzt:nzt+1,14,tn) = sums_l(nzt:nzt+1,14,tn) + &
827                                    vswst(j,i) * rmask(j,i,sr)    ! w"v"
828                sums_l(nzt:nzt+1,16,tn) = sums_l(nzt:nzt+1,16,tn) + &
829                                    tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
830                sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + &
831                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
832                sums_l(nzt:nzt+1,61,tn) = sums_l(nzt:nzt+1,61,tn) + &
833                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
834
835                IF ( ocean )  THEN
836                   sums_l(nzt,65,tn) = sums_l(nzt,65,tn) + &
837                                       saswst(j,i) * rmask(j,i,sr)  ! w"sa"
838                ENDIF
839                IF ( humidity )  THEN
840                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) +                     &
841                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
842                   sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                   &
843                                       ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) *     &
844                                       tswst(j,i) + 0.61_wp * pt(nzt,j,i) *    &
845                                                             qswst(j,i) )
846                   IF ( cloud_droplets )  THEN
847                      sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                &
848                                          ( 1.0_wp + 0.61_wp * q(nzt,j,i) -    &
849                                            ql(nzt,j,i) ) * tswst(j,i) +       &
850                                            0.61_wp * pt(nzt,j,i) * qswst(j,i) )
851                   ENDIF
852                   IF ( cloud_physics )  THEN
853!
854!--                   Formula does not work if ql(nzb) /= 0.0
855                      sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + &   ! w"q" (w"qv")
856                                          qswst(j,i) * rmask(j,i,sr)
857                   ENDIF
858                ENDIF
859                IF ( passive_scalar )  THEN
860                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
861                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
862                ENDIF
863             ENDIF
864
865!
866!--          Resolved fluxes (can be computed for all horizontal points)
867!--          NOTE: for simplicity, nzb_s_inner is used below, although strictly
868!--          ----  speaking the following k-loop would have to be split up and
869!--                rearranged according to the staggered grid.
870             DO  k = nzb_s_inner(j,i), nzt
871                ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +                  &
872                                 u(k+1,j,i) - hom(k+1,1,1,sr) )
873                vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +                  &
874                                 v(k+1,j,i) - hom(k+1,1,2,sr) )
875                pts = 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) +                 &
876                                 pt(k+1,j,i) - hom(k+1,1,4,sr) )
877
878!--             Higher moments
879                sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 *        &
880                                                    rmask(j,i,sr)
881                sums_l(k,36,tn) = sums_l(k,36,tn) + pts**2 * w(k,j,i) *        &
882                                                    rmask(j,i,sr)
883
884!
885!--             Salinity flux and density (density does not belong to here,
886!--             but so far there is no other suitable place to calculate)
887                IF ( ocean )  THEN
888                   IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
889                      pts = 0.5_wp * ( sa(k,j,i)   - hom(k,1,23,sr) +          &
890                                       sa(k+1,j,i) - hom(k+1,1,23,sr) )
891                      sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) *     &
892                                        rmask(j,i,sr)
893                   ENDIF
894                   sums_l(k,64,tn) = sums_l(k,64,tn) + rho(k,j,i) *            &
895                                                       rmask(j,i,sr)
896                   sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) *           &
897                                                       rmask(j,i,sr)
898                ENDIF
899
900!
901!--             Buoyancy flux, water flux, humidity flux, liquid water
902!--             content, rain drop concentration and rain water content
903                IF ( humidity )  THEN
904                   IF ( cloud_physics .OR. cloud_droplets )  THEN
905                      pts = 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +         &
906                                    vpt(k+1,j,i) - hom(k+1,1,44,sr) )
907                      sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *     &
908                                                          rmask(j,i,sr)
909                      sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * rmask(j,i,sr)
910
911                      IF ( .NOT. cloud_droplets )  THEN
912                         pts = 0.5_wp *                                        &
913                              ( ( q(k,j,i) - ql(k,j,i) ) -                     &
914                              hom(k,1,42,sr) +                                 &
915                              ( q(k+1,j,i) - ql(k+1,j,i) ) -                   &
916                              hom(k+1,1,42,sr) )
917                         sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) *  &
918                                                             rmask(j,i,sr)
919                         sums_l(k,75,tn) = sums_l(k,75,tn) + qc(k,j,i) *       &
920                                                             rmask(j,i,sr)
921                         sums_l(k,76,tn) = sums_l(k,76,tn) + prr(k,j,i) *      &
922                                                             rmask(j,i,sr)
923                         IF ( microphysics_seifert )  THEN
924                            sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) *    &
925                                                                rmask(j,i,sr)
926                            sums_l(k,74,tn) = sums_l(k,74,tn) + qr(k,j,i) *    &
927                                                                rmask(j,i,sr)
928                         ENDIF
929                      ENDIF
930
931                   ELSE
932                      IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
933                         pts = 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +      &
934                                          vpt(k+1,j,i) - hom(k+1,1,44,sr) )
935                         sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *  &
936                                                             rmask(j,i,sr)
937                      ELSE IF ( ws_scheme_sca .AND. sr == 0 )  THEN
938                         sums_l(k,46,tn) = ( 1.0_wp + 0.61_wp *                & 
939                                             hom(k,1,41,sr) ) *                &
940                                           sums_l(k,17,tn) +                   &
941                                           0.61_wp * hom(k,1,4,sr) *           &
942                                           sums_l(k,49,tn)
943                      END IF
944                   END IF
945                ENDIF
946!
947!--             Passive scalar flux
948                IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca                &
949                     .OR. sr /= 0 ) )  THEN
950                   pts = 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +              &
951                                    q(k+1,j,i) - hom(k+1,1,41,sr) )
952                   sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *        &
953                                                       rmask(j,i,sr)
954                ENDIF
955
956!
957!--             Energy flux w*e*
958!--             has to be adjusted
959                sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5_wp *        &
960                                             ( ust**2 + vst**2 + w(k,j,i)**2 ) &
961                                             * rmask(j,i,sr)
962             ENDDO
963          ENDDO
964       ENDDO
965!
966!--    For speed optimization fluxes which have been computed in part directly
967!--    inside the WS advection routines are treated seperatly
968!--    Momentum fluxes first:
969       IF ( .NOT. ws_scheme_mom .OR. sr /= 0  )  THEN
970         !$OMP DO
971         DO  i = nxl, nxr
972            DO  j = nys, nyn
973               DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
974                  ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +                &
975                                   u(k+1,j,i) - hom(k+1,1,1,sr) )
976                  vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +                &
977                                   v(k+1,j,i) - hom(k+1,1,2,sr) )
978!
979!--               Momentum flux w*u*
980                  sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5_wp *                 &
981                                                    ( w(k,j,i-1) + w(k,j,i) )  &
982                                                    * ust * rmask(j,i,sr)
983!
984!--               Momentum flux w*v*
985                  sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5_wp *                 &
986                                                    ( w(k,j-1,i) + w(k,j,i) )  &
987                                                    * vst * rmask(j,i,sr)
988               ENDDO
989            ENDDO
990         ENDDO
991
992       ENDIF
993       IF ( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
994         !$OMP DO
995         DO  i = nxl, nxr
996            DO  j = nys, nyn
997               DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
998!
999!--               Vertical heat flux
1000                  sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5_wp *                 &
1001                           ( pt(k,j,i)   - hom(k,1,4,sr) +                     &
1002                             pt(k+1,j,i) - hom(k+1,1,4,sr) )                   &
1003                           * w(k,j,i) * rmask(j,i,sr)
1004                  IF ( humidity )  THEN
1005                     pts = 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +            &
1006                                      q(k+1,j,i) - hom(k+1,1,41,sr) )
1007                     sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *      &
1008                                       rmask(j,i,sr)
1009                  ENDIF
1010               ENDDO
1011            ENDDO
1012         ENDDO
1013
1014       ENDIF
1015
1016!
1017!--    Density at top follows Neumann condition
1018       IF ( ocean )  THEN
1019          sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
1020          sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
1021       ENDIF
1022
1023!
1024!--    Divergence of vertical flux of resolved scale energy and pressure
1025!--    fluctuations as well as flux of pressure fluctuation itself (68).
1026!--    First calculate the products, then the divergence.
1027!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
1028       IF ( hom(nzb+1,2,55,0) /= 0.0_wp  .OR.  hom(nzb+1,2,68,0) /= 0.0_wp )   &
1029       THEN
1030          sums_ll = 0.0_wp  ! local array
1031
1032          !$OMP DO
1033          DO  i = nxl, nxr
1034             DO  j = nys, nyn
1035                DO  k = nzb_s_inner(j,i)+1, nzt
1036
1037                   sums_ll(k,1) = sums_ll(k,1) + 0.5_wp * w(k,j,i) * (         &
1038                  ( 0.25_wp * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1) )  &
1039                            - 0.5_wp * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) ) )**2&
1040                + ( 0.25_wp * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i) )  &
1041                            - 0.5_wp * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) ) )**2&
1042                + w(k,j,i)**2                                        )
1043
1044                   sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i)             &
1045                                               * ( p(k,j,i) + p(k+1,j,i) )
1046
1047                ENDDO
1048             ENDDO
1049          ENDDO
1050          sums_ll(0,1)     = 0.0_wp    ! because w is zero at the bottom
1051          sums_ll(nzt+1,1) = 0.0_wp
1052          sums_ll(0,2)     = 0.0_wp
1053          sums_ll(nzt+1,2) = 0.0_wp
1054
1055          DO  k = nzb+1, nzt
1056             sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
1057             sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
1058             sums_l(k,68,tn) = sums_ll(k,2)
1059          ENDDO
1060          sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
1061          sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
1062          sums_l(nzb,68,tn) = 0.0_wp    ! because w* = 0 at nzb
1063
1064       ENDIF
1065
1066!
1067!--    Divergence of vertical flux of SGS TKE and the flux itself (69)
1068       IF ( hom(nzb+1,2,57,0) /= 0.0_wp  .OR.  hom(nzb+1,2,69,0) /= 0.0_wp )   &
1069       THEN
1070          !$OMP DO
1071          DO  i = nxl, nxr
1072             DO  j = nys, nyn
1073                DO  k = nzb_s_inner(j,i)+1, nzt
1074
1075                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5_wp * (              &
1076                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
1077                 - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
1078                                                                ) * ddzw(k)
1079
1080                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5_wp * (              &
1081                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
1082                                                                )
1083
1084                ENDDO
1085             ENDDO
1086          ENDDO
1087          sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
1088          sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
1089
1090       ENDIF
1091
1092!
1093!--    Horizontal heat fluxes (subgrid, resolved, total).
1094!--    Do it only, if profiles shall be plotted.
1095       IF ( hom(nzb+1,2,58,0) /= 0.0_wp ) THEN
1096
1097          !$OMP DO
1098          DO  i = nxl, nxr
1099             DO  j = nys, nyn
1100                DO  k = nzb_s_inner(j,i)+1, nzt
1101!
1102!--                Subgrid horizontal heat fluxes u"pt", v"pt"
1103                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5_wp *                &
1104                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
1105                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
1106                                                 * ddx * rmask(j,i,sr)
1107                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
1108                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
1109                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
1110                                                 * ddy * rmask(j,i,sr)
1111!
1112!--                Resolved horizontal heat fluxes u*pt*, v*pt*
1113                   sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
1114                                                  ( u(k,j,i) - hom(k,1,1,sr) ) &
1115                                    * 0.5_wp * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
1116                                                 pt(k,j,i)   - hom(k,1,4,sr) )
1117                   pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
1118                                    pt(k,j,i)   - hom(k,1,4,sr) )
1119                   sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
1120                                                  ( v(k,j,i) - hom(k,1,2,sr) ) &
1121                                    * 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
1122                                                 pt(k,j,i)   - hom(k,1,4,sr) )
1123                ENDDO
1124             ENDDO
1125          ENDDO
1126!
1127!--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
1128          sums_l(nzb,58,tn) = 0.0_wp
1129          sums_l(nzb,59,tn) = 0.0_wp
1130          sums_l(nzb,60,tn) = 0.0_wp
1131          sums_l(nzb,61,tn) = 0.0_wp
1132          sums_l(nzb,62,tn) = 0.0_wp
1133          sums_l(nzb,63,tn) = 0.0_wp
1134
1135       ENDIF
1136
1137!
1138!--    Collect current large scale advection and subsidence tendencies for
1139!--    data output
1140       IF ( large_scale_forcing  .AND.  ( simulated_time > 0.0_wp ) )  THEN
1141!
1142!--       Interpolation in time of LSF_DATA
1143          nt = 1
1144          DO WHILE ( simulated_time - dt_3d > time_vert(nt) )
1145             nt = nt + 1
1146          ENDDO
1147          IF ( simulated_time - dt_3d /= time_vert(nt) )  THEN
1148            nt = nt - 1
1149          ENDIF
1150
1151          fac = ( simulated_time - dt_3d - time_vert(nt) )                     &
1152                / ( time_vert(nt+1)-time_vert(nt) )
1153
1154
1155          DO  k = nzb, nzt
1156             sums_ls_l(k,0) = td_lsa_lpt(k,nt)                                 &
1157                              + fac * ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )
1158             sums_ls_l(k,1) = td_lsa_q(k,nt)                                   &
1159                              + fac * ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) )
1160          ENDDO
1161
1162          sums_ls_l(nzt+1,0) = sums_ls_l(nzt,0)
1163          sums_ls_l(nzt+1,1) = sums_ls_l(nzt,1)
1164
1165          IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
1166
1167             DO  k = nzb, nzt
1168                sums_ls_l(k,2) = td_sub_lpt(k,nt) + fac *                      &
1169                                 ( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )
1170                sums_ls_l(k,3) = td_sub_q(k,nt) + fac *                        &
1171                                 ( td_sub_q(k,nt+1) - td_sub_q(k,nt) )
1172             ENDDO
1173
1174             sums_ls_l(nzt+1,2) = sums_ls_l(nzt,2)
1175             sums_ls_l(nzt+1,3) = sums_ls_l(nzt,3)
1176
1177          ENDIF
1178
1179       ENDIF
1180
1181
1182       IF ( land_surface )  THEN
1183          !$OMP DO
1184          DO  i = nxl, nxr
1185             DO  j =  nys, nyn
1186                DO  k = nzb_soil, nzt_soil
1187                   sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil(k,j,i)         &
1188                                      * rmask(j,i,sr)
1189                   sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil(k,j,i)         &
1190                                      * rmask(j,i,sr)
1191                ENDDO
1192             ENDDO
1193          ENDDO
1194       ENDIF
1195       
1196       IF ( radiation .AND. radiation_scheme == 'rrtmg' )  THEN
1197          !$OMP DO
1198          DO  i = nxl, nxr
1199             DO  j =  nys, nyn
1200                DO  k = nzb_s_inner(j,i)+1, nzt+1
1201                   sums_l(k,102,tn)  = sums_l(k,102,tn)  + rad_lw_in(k,j,i)    &
1202                                       * rmask(j,i,sr)
1203                   sums_l(k,103,tn)  = sums_l(k,103,tn)  + rad_lw_out(k,j,i)   &
1204                                       * rmask(j,i,sr)
1205                   sums_l(k,104,tn)  = sums_l(k,104,tn)  + rad_sw_in(k,j,i)    &
1206                                       * rmask(j,i,sr)
1207                   sums_l(k,105,tn)  = sums_l(k,105,tn)  + rad_sw_out(k,j,i)   &
1208                                       * rmask(j,i,sr)
1209                   sums_l(k,106,tn)  = sums_l(k,106,tn)  + rad_lw_cs_hr(k,j,i) &
1210                                       * rmask(j,i,sr)
1211                   sums_l(k,107,tn)  = sums_l(k,107,tn)  + rad_lw_hr(k,j,i)    &
1212                                       * rmask(j,i,sr)
1213                   sums_l(k,108,tn)  = sums_l(k,108,tn)  + rad_sw_cs_hr(k,j,i) &
1214                                       * rmask(j,i,sr)
1215                   sums_l(k,109,tn)  = sums_l(k,109,tn)  + rad_sw_hr(k,j,i)    &
1216                                       * rmask(j,i,sr)
1217                ENDDO
1218             ENDDO
1219          ENDDO
1220       ENDIF
1221!
1222!--    Calculate the user-defined profiles
1223       CALL user_statistics( 'profiles', sr, tn )
1224       !$OMP END PARALLEL
1225
1226!
1227!--    Summation of thread sums
1228       IF ( threads_per_task > 1 )  THEN
1229          DO  i = 1, threads_per_task-1
1230             sums_l(:,3,0)          = sums_l(:,3,0) + sums_l(:,3,i)
1231             sums_l(:,4:40,0)       = sums_l(:,4:40,0) + sums_l(:,4:40,i)
1232             sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
1233                                      sums_l(:,45:pr_palm,i)
1234             IF ( max_pr_user > 0 )  THEN
1235                sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
1236                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
1237                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
1238             ENDIF
1239          ENDDO
1240       ENDIF
1241
1242#if defined( __parallel )
1243
1244!
1245!--    Compute total sum from local sums
1246       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1247       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL,   &
1248                           MPI_SUM, comm2d, ierr )
1249       IF ( large_scale_forcing )  THEN
1250          CALL MPI_ALLREDUCE( sums_ls_l(nzb,2), sums(nzb,83), ngp_sums_ls,     &
1251                              MPI_REAL, MPI_SUM, comm2d, ierr )
1252       ENDIF
1253#else
1254       sums = sums_l(:,:,0)
1255       IF ( large_scale_forcing )  THEN
1256          sums(:,81:88) = sums_ls_l
1257       ENDIF
1258#endif
1259
1260!
1261!--    Final values are obtained by division by the total number of grid points
1262!--    used for summation. After that store profiles.
1263!--    Check, if statistical regions do contain at least one grid point at the
1264!--    respective k-level, otherwise division by zero will lead to undefined
1265!--    values, which may cause e.g. problems with NetCDF output
1266!--    Profiles:
1267       DO  k = nzb, nzt+1
1268          sums(k,3)             = sums(k,3)             / ngp_2dh(sr)
1269          sums(k,12:22)         = sums(k,12:22)         / ngp_2dh(sr)
1270          sums(k,30:32)         = sums(k,30:32)         / ngp_2dh(sr)
1271          sums(k,35:39)         = sums(k,35:39)         / ngp_2dh(sr)
1272          sums(k,45:53)         = sums(k,45:53)         / ngp_2dh(sr)
1273          sums(k,55:63)         = sums(k,55:63)         / ngp_2dh(sr)
1274          sums(k,81:88)         = sums(k,81:88)         / ngp_2dh(sr)
1275          sums(k,89:114)        = sums(k,89:114)        / ngp_2dh(sr)
1276          IF ( ngp_2dh_s_inner(k,sr) /= 0 )  THEN
1277             sums(k,8:11)          = sums(k,8:11)          / ngp_2dh_s_inner(k,sr)
1278             sums(k,23:29)         = sums(k,23:29)         / ngp_2dh_s_inner(k,sr)
1279             sums(k,33:34)         = sums(k,33:34)         / ngp_2dh_s_inner(k,sr)
1280             sums(k,40)            = sums(k,40)            / ngp_2dh_s_inner(k,sr)
1281             sums(k,54)            = sums(k,54)            / ngp_2dh_s_inner(k,sr)
1282             sums(k,64)            = sums(k,64)            / ngp_2dh_s_inner(k,sr)
1283             sums(k,70:80)         = sums(k,70:80)         / ngp_2dh_s_inner(k,sr)
1284             sums(k,115:pr_palm-2) = sums(k,115:pr_palm-2) / ngp_2dh_s_inner(k,sr)
1285          ENDIF
1286       ENDDO
1287
1288!--    u* and so on
1289!--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
1290!--    size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
1291!--    above the topography, they are being divided by ngp_2dh(sr)
1292       sums(nzb:nzb+3,pr_palm)    = sums(nzb:nzb+3,pr_palm)    / &
1293                                    ngp_2dh(sr)
1294       sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)       / &    ! qs
1295                                    ngp_2dh(sr)
1296!--    eges, e*
1297       sums(nzb+4:nzb+5,pr_palm)  = sums(nzb+4:nzb+5,pr_palm)  / &
1298                                    ngp_3d(sr)
1299!--    Old and new divergence
1300       sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
1301                                    ngp_3d_inner(sr)
1302
1303!--    User-defined profiles
1304       IF ( max_pr_user > 0 )  THEN
1305          DO  k = nzb, nzt+1
1306             sums(k,pr_palm+1:pr_palm+max_pr_user) = &
1307                                    sums(k,pr_palm+1:pr_palm+max_pr_user) / &
1308                                    ngp_2dh_s_inner(k,sr)
1309          ENDDO
1310       ENDIF
1311
1312!
1313!--    Collect horizontal average in hom.
1314!--    Compute deduced averages (e.g. total heat flux)
1315       hom(:,1,3,sr)  = sums(:,3)      ! w
1316       hom(:,1,8,sr)  = sums(:,8)      ! e     profiles 5-7 are initial profiles
1317       hom(:,1,9,sr)  = sums(:,9)      ! km
1318       hom(:,1,10,sr) = sums(:,10)     ! kh
1319       hom(:,1,11,sr) = sums(:,11)     ! l
1320       hom(:,1,12,sr) = sums(:,12)     ! w"u"
1321       hom(:,1,13,sr) = sums(:,13)     ! w*u*
1322       hom(:,1,14,sr) = sums(:,14)     ! w"v"
1323       hom(:,1,15,sr) = sums(:,15)     ! w*v*
1324       hom(:,1,16,sr) = sums(:,16)     ! w"pt"
1325       hom(:,1,17,sr) = sums(:,17)     ! w*pt*
1326       hom(:,1,18,sr) = sums(:,16) + sums(:,17)    ! wpt
1327       hom(:,1,19,sr) = sums(:,12) + sums(:,13)    ! wu
1328       hom(:,1,20,sr) = sums(:,14) + sums(:,15)    ! wv
1329       hom(:,1,21,sr) = sums(:,21)     ! w*pt*BC
1330       hom(:,1,22,sr) = sums(:,16) + sums(:,21)    ! wptBC
1331                                       ! profile 24 is initial profile (sa)
1332                                       ! profiles 25-29 left empty for initial
1333                                       ! profiles
1334       hom(:,1,30,sr) = sums(:,30)     ! u*2
1335       hom(:,1,31,sr) = sums(:,31)     ! v*2
1336       hom(:,1,32,sr) = sums(:,32)     ! w*2
1337       hom(:,1,33,sr) = sums(:,33)     ! pt*2
1338       hom(:,1,34,sr) = sums(:,34)     ! e*
1339       hom(:,1,35,sr) = sums(:,35)     ! w*2pt*
1340       hom(:,1,36,sr) = sums(:,36)     ! w*pt*2
1341       hom(:,1,37,sr) = sums(:,37)     ! w*e*
1342       hom(:,1,38,sr) = sums(:,38)     ! w*3
1343       hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20_wp )**1.5_wp   ! Sw
1344       hom(:,1,40,sr) = sums(:,40)     ! p
1345       hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
1346       hom(:,1,46,sr) = sums(:,46)     ! w*vpt*       
1347       hom(:,1,47,sr) = sums(:,45) + sums(:,46)    ! wvpt
1348       hom(:,1,48,sr) = sums(:,48)     ! w"q" (w"qv")
1349       hom(:,1,49,sr) = sums(:,49)     ! w*q* (w*qv*)
1350       hom(:,1,50,sr) = sums(:,48) + sums(:,49)    ! wq (wqv)
1351       hom(:,1,51,sr) = sums(:,51)     ! w"qv"
1352       hom(:,1,52,sr) = sums(:,52)     ! w*qv*       
1353       hom(:,1,53,sr) = sums(:,52) + sums(:,51)    ! wq (wqv)
1354       hom(:,1,54,sr) = sums(:,54)     ! ql
1355       hom(:,1,55,sr) = sums(:,55)     ! w*u*u*/dz
1356       hom(:,1,56,sr) = sums(:,56)     ! w*p*/dz
1357       hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho )/dz
1358       hom(:,1,58,sr) = sums(:,58)     ! u"pt"
1359       hom(:,1,59,sr) = sums(:,59)     ! u*pt*
1360       hom(:,1,60,sr) = sums(:,58) + sums(:,59)    ! upt_t
1361       hom(:,1,61,sr) = sums(:,61)     ! v"pt"
1362       hom(:,1,62,sr) = sums(:,62)     ! v*pt*
1363       hom(:,1,63,sr) = sums(:,61) + sums(:,62)    ! vpt_t
1364       hom(:,1,64,sr) = sums(:,64)     ! rho
1365       hom(:,1,65,sr) = sums(:,65)     ! w"sa"
1366       hom(:,1,66,sr) = sums(:,66)     ! w*sa*
1367       hom(:,1,67,sr) = sums(:,65) + sums(:,66)    ! wsa
1368       hom(:,1,68,sr) = sums(:,68)     ! w*p*
1369       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho
1370       hom(:,1,70,sr) = sums(:,70)     ! q*2
1371       hom(:,1,71,sr) = sums(:,71)     ! prho
1372       hom(:,1,72,sr) = hyp * 1E-4_wp  ! hyp in dbar
1373       hom(:,1,73,sr) = sums(:,73)     ! nr
1374       hom(:,1,74,sr) = sums(:,74)     ! qr
1375       hom(:,1,75,sr) = sums(:,75)     ! qc
1376       hom(:,1,76,sr) = sums(:,76)     ! prr (precipitation rate)
1377                                       ! 77 is initial density profile
1378       hom(:,1,78,sr) = ug             ! ug
1379       hom(:,1,79,sr) = vg             ! vg
1380       hom(:,1,80,sr) = w_subs         ! w_subs
1381
1382       IF ( large_scale_forcing )  THEN
1383          hom(:,1,81,sr) = sums_ls_l(:,0)          ! td_lsa_lpt
1384          hom(:,1,82,sr) = sums_ls_l(:,1)          ! td_lsa_q
1385          IF ( use_subsidence_tendencies )  THEN
1386             hom(:,1,83,sr) = sums_ls_l(:,2)       ! td_sub_lpt
1387             hom(:,1,84,sr) = sums_ls_l(:,3)       ! td_sub_q
1388          ELSE
1389             hom(:,1,83,sr) = sums(:,83)           ! td_sub_lpt
1390             hom(:,1,84,sr) = sums(:,84)           ! td_sub_q
1391          ENDIF
1392          hom(:,1,85,sr) = sums(:,85)              ! td_nud_lpt
1393          hom(:,1,86,sr) = sums(:,86)              ! td_nud_q
1394          hom(:,1,87,sr) = sums(:,87)              ! td_nud_u
1395          hom(:,1,88,sr) = sums(:,88)              ! td_nud_v
1396       ENDIF
1397
1398       IF ( land_surface )  THEN
1399          hom(:,1,89,sr) = sums(:,89)              ! t_soil
1400                                                   ! 90 is initial t_soil profile
1401          hom(:,1,91,sr) = sums(:,91)              ! m_soil
1402                                                   ! 92 is initial m_soil profile
1403          hom(:,1,93,sr)  = sums(:,93)             ! ghf_eb
1404          hom(:,1,94,sr)  = sums(:,94)             ! shf_eb
1405          hom(:,1,95,sr)  = sums(:,95)             ! qsws_eb
1406          hom(:,1,96,sr)  = sums(:,96)             ! qsws_liq_eb
1407          hom(:,1,97,sr)  = sums(:,97)             ! qsws_soil_eb
1408          hom(:,1,98,sr)  = sums(:,98)             ! qsws_veg_eb
1409          hom(:,1,99,sr)  = sums(:,99)             ! r_a
1410          hom(:,1,100,sr) = sums(:,100)            ! r_s
1411
1412       ENDIF
1413
1414       IF ( radiation )  THEN
1415          hom(:,1,101,sr) = sums(:,101)            ! rad_net
1416          hom(:,1,102,sr) = sums(:,102)            ! rad_lw_in
1417          hom(:,1,103,sr) = sums(:,103)            ! rad_lw_out
1418          hom(:,1,104,sr) = sums(:,104)            ! rad_sw_in
1419          hom(:,1,105,sr) = sums(:,105)            ! rad_sw_out
1420
1421          IF ( radiation_scheme == 'rrtmg' )  THEN
1422#if defined ( __rrtmg )
1423             hom(:,1,106,sr) = sums(:,106)            ! rad_lw_cs_hr
1424             hom(:,1,107,sr) = sums(:,107)            ! rad_lw_hr
1425             hom(:,1,108,sr) = sums(:,108)            ! rad_sw_cs_hr
1426             hom(:,1,109,sr) = sums(:,109)            ! rad_sw_hr
1427
1428             hom(:,1,110,sr) = sums(:,110)            ! rrtm_aldif
1429             hom(:,1,111,sr) = sums(:,111)            ! rrtm_aldir
1430             hom(:,1,112,sr) = sums(:,112)            ! rrtm_asdif
1431             hom(:,1,113,sr) = sums(:,113)            ! rrtm_asdir
1432#endif
1433          ENDIF
1434
1435
1436       ENDIF
1437
1438       hom(:,1,114,sr) = sums(:,114)            !: L
1439
1440       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
1441                                       ! u*, w'u', w'v', t* (in last profile)
1442
1443       IF ( max_pr_user > 0 )  THEN    ! user-defined profiles
1444          hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
1445                               sums(:,pr_palm+1:pr_palm+max_pr_user)
1446       ENDIF
1447
1448!
1449!--    Determine the boundary layer height using two different schemes.
1450!--    First scheme: Starting from the Earth's (Ocean's) surface, look for the
1451!--    first relative minimum (maximum) of the total heat flux.
1452!--    The corresponding height is assumed as the boundary layer height, if it
1453!--    is less than 1.5 times the height where the heat flux becomes negative
1454!--    (positive) for the first time.
1455       z_i(1) = 0.0_wp
1456       first = .TRUE.
1457
1458       IF ( ocean )  THEN
1459          DO  k = nzt, nzb+1, -1
1460             IF ( first  .AND.  hom(k,1,18,sr) < -1.0E-8_wp )  THEN
1461                first = .FALSE.
1462                height = zw(k)
1463             ENDIF
1464             IF ( hom(k,1,18,sr) < -1.0E-8_wp  .AND.                           &
1465                  hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
1466                IF ( zw(k) < 1.5_wp * height )  THEN
1467                   z_i(1) = zw(k)
1468                ELSE
1469                   z_i(1) = height
1470                ENDIF
1471                EXIT
1472             ENDIF
1473          ENDDO
1474       ELSE
1475          DO  k = nzb, nzt-1
1476             IF ( first  .AND.  hom(k,1,18,sr) < -1.0E-8_wp )  THEN
1477                first = .FALSE.
1478                height = zw(k)
1479             ENDIF
1480             IF ( hom(k,1,18,sr) < -1.0E-8_wp  .AND.                           &
1481                  hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
1482                IF ( zw(k) < 1.5_wp * height )  THEN
1483                   z_i(1) = zw(k)
1484                ELSE
1485                   z_i(1) = height
1486                ENDIF
1487                EXIT
1488             ENDIF
1489          ENDDO
1490       ENDIF
1491
1492!
1493!--    Second scheme: Gradient scheme from Sullivan et al. (1998), modified
1494!--    by Uhlenbrock(2006). The boundary layer height is the height with the
1495!--    maximal local temperature gradient: starting from the second (the last
1496!--    but one) vertical gridpoint, the local gradient must be at least
1497!--    0.2K/100m and greater than the next four gradients.
1498!--    WARNING: The threshold value of 0.2K/100m must be adjusted for the
1499!--             ocean case!
1500       z_i(2) = 0.0_wp
1501       DO  k = nzb+1, nzt+1
1502          dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
1503       ENDDO
1504       dptdz_threshold = 0.2_wp / 100.0_wp
1505
1506       IF ( ocean )  THEN
1507          DO  k = nzt+1, nzb+5, -1
1508             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
1509                  dptdz(k) > dptdz(k-1)  .AND.  dptdz(k) > dptdz(k-2)  .AND.  &
1510                  dptdz(k) > dptdz(k-3)  .AND.  dptdz(k) > dptdz(k-4) )  THEN
1511                z_i(2) = zw(k-1)
1512                EXIT
1513             ENDIF
1514          ENDDO
1515       ELSE
1516          DO  k = nzb+1, nzt-3
1517             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
1518                  dptdz(k) > dptdz(k+1)  .AND.  dptdz(k) > dptdz(k+2)  .AND.  &
1519                  dptdz(k) > dptdz(k+3)  .AND.  dptdz(k) > dptdz(k+4) )  THEN
1520                z_i(2) = zw(k-1)
1521                EXIT
1522             ENDIF
1523          ENDDO
1524       ENDIF
1525
1526       hom(nzb+6,1,pr_palm,sr) = z_i(1)
1527       hom(nzb+7,1,pr_palm,sr) = z_i(2)
1528
1529!
1530!--    Determine vertical index which is nearest to the mean surface level
1531!--    height of the respective statistic region
1532       DO  k = nzb, nzt
1533          IF ( zw(k) >= mean_surface_level_height(sr) )  THEN
1534             k_surface_level = k
1535             EXIT
1536          ENDIF
1537       ENDDO
1538!
1539!--    Computation of both the characteristic vertical velocity and
1540!--    the characteristic convective boundary layer temperature.
1541!--    The inversion height entering into the equation is defined with respect
1542!--    to the mean surface level height of the respective statistic region.
1543!--    The horizontal average at surface level index + 1 is input for the
1544!--    average temperature.
1545       IF ( hom(k_surface_level,1,18,sr) > 1.0E-8_wp  .AND.  z_i(1) /= 0.0_wp )&
1546       THEN
1547          hom(nzb+8,1,pr_palm,sr) = &
1548             ( g / hom(k_surface_level+1,1,4,sr) * hom(k_surface_level,1,18,sr)&
1549             * ABS( z_i(1) - mean_surface_level_height(sr) ) )**0.333333333_wp
1550       ELSE
1551          hom(nzb+8,1,pr_palm,sr)  = 0.0_wp
1552       ENDIF
1553
1554!
1555!--    Collect the time series quantities
1556       ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)     ! E
1557       ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)     ! E*
1558       ts_value(3,sr) = dt_3d
1559       ts_value(4,sr) = hom(nzb,1,pr_palm,sr)       ! u*
1560       ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr)     ! th*
1561       ts_value(6,sr) = u_max
1562       ts_value(7,sr) = v_max
1563       ts_value(8,sr) = w_max
1564       ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr)    ! new divergence
1565       ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr)    ! old Divergence
1566       ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr)    ! z_i(1)
1567       ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr)    ! z_i(2)
1568       ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr)    ! w*
1569       ts_value(14,sr) = hom(nzb,1,16,sr)           ! w'pt'   at k=0
1570       ts_value(15,sr) = hom(nzb+1,1,16,sr)         ! w'pt'   at k=1
1571       ts_value(16,sr) = hom(nzb+1,1,18,sr)         ! wpt     at k=1
1572       ts_value(17,sr) = hom(nzb,1,4,sr)            ! pt(0)
1573       ts_value(18,sr) = hom(nzb+1,1,4,sr)          ! pt(zp)
1574       ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr)    ! u'w'    at k=0
1575       ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)    ! v'w'    at k=0
1576       ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
1577
1578       IF ( .NOT. neutral )  THEN
1579          ts_value(22,sr) = hom(nzb,1,114,sr)          ! L
1580       ELSE
1581          ts_value(22,sr) = 1.0E10_wp
1582       ENDIF
1583
1584       ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
1585
1586!
1587!--    Collect land surface model timeseries
1588       IF ( land_surface )  THEN
1589          ts_value(dots_soil  ,sr) = hom(nzb,1,93,sr)           ! ghf_eb
1590          ts_value(dots_soil+1,sr) = hom(nzb,1,94,sr)           ! shf_eb
1591          ts_value(dots_soil+2,sr) = hom(nzb,1,95,sr)           ! qsws_eb
1592          ts_value(dots_soil+3,sr) = hom(nzb,1,96,sr)           ! qsws_liq_eb
1593          ts_value(dots_soil+4,sr) = hom(nzb,1,97,sr)           ! qsws_soil_eb
1594          ts_value(dots_soil+5,sr) = hom(nzb,1,98,sr)           ! qsws_veg_eb
1595          ts_value(dots_soil+6,sr) = hom(nzb,1,99,sr)           ! r_a
1596          ts_value(dots_soil+7,sr) = hom(nzb,1,100,sr)          ! r_s
1597       ENDIF
1598!
1599!--    Collect radiation model timeseries
1600       IF ( radiation )  THEN
1601          ts_value(dots_rad,sr)   = hom(nzb,1,101,sr)          ! rad_net
1602          ts_value(dots_rad+1,sr) = hom(nzb,1,102,sr)          ! rad_lw_in
1603          ts_value(dots_rad+2,sr) = hom(nzb,1,103,sr)          ! rad_lw_out
1604          ts_value(dots_rad+3,sr) = hom(nzb,1,104,sr)          ! rad_sw_in
1605          ts_value(dots_rad+4,sr) = hom(nzb,1,105,sr)          ! rad_sw_out
1606
1607#if defined ( __rrtmg )
1608          IF ( radiation_scheme == 'rrtmg' )  THEN
1609             ts_value(dots_rad+5,sr) = hom(nzb,1,110,sr)          ! rrtm_aldif
1610             ts_value(dots_rad+6,sr) = hom(nzb,1,111,sr)          ! rrtm_aldir
1611             ts_value(dots_rad+7,sr) = hom(nzb,1,112,sr)          ! rrtm_asdif
1612             ts_value(dots_rad+8,sr) = hom(nzb,1,113,sr)          ! rrtm_asdir
1613          ENDIF
1614#endif
1615
1616       ENDIF
1617
1618!
1619!--    Calculate additional statistics provided by the user interface
1620       CALL user_statistics( 'time_series', sr, 0 )
1621
1622    ENDDO    ! loop of the subregions
1623
1624!
1625!-- If required, sum up horizontal averages for subsequent time averaging
1626    IF ( do_sum )  THEN
1627       IF ( average_count_pr == 0 )  hom_sum = 0.0_wp
1628       hom_sum = hom_sum + hom(:,1,:,:)
1629       average_count_pr = average_count_pr + 1
1630       do_sum = .FALSE.
1631    ENDIF
1632
1633!
1634!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
1635!-- This flag is reset after each time step in time_integration.
1636    flow_statistics_called = .TRUE.
1637
1638    CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
1639
1640
1641 END SUBROUTINE flow_statistics
1642
1643
1644#else
1645
1646
1647!------------------------------------------------------------------------------!
1648! Description:
1649! ------------
1650!> flow statistics - accelerator version
1651!------------------------------------------------------------------------------!
1652 SUBROUTINE flow_statistics
1653
1654    USE arrays_3d,                                                             &
1655        ONLY:  ddzu, ddzw, e, hyp, km, kh, nr, p, prho, pt, q, qc, ql, qr, qs, &
1656               qsws, qswst, rho, sa, saswsb, saswst, shf, td_lsa_lpt, td_lsa_q,&
1657               td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us, usws,    &
1658               uswst, vsws, v, vg, vpt, vswst, w, w_subs, zw
1659                 
1660       
1661    USE cloud_parameters,                                                      &
1662        ONLY:  l_d_cp, prr, pt_d_t
1663       
1664    USE control_parameters,                                                    &
1665        ONLY :  average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
1666                dt_3d, g, humidity, kappa, large_scale_forcing,                &
1667                large_scale_subsidence, max_pr_user, message_string,           &
1668                microphysics_seifert, neutral, ocean, passive_scalar,          &
1669                simulated_time, use_subsidence_tendencies, use_surface_fluxes, &
1670                use_top_fluxes, ws_scheme_mom, ws_scheme_sca
1671       
1672    USE cpulog,                                                                &
1673        ONLY:  cpu_log, log_point
1674       
1675    USE grid_variables,                                                        &
1676        ONLY:  ddx, ddy
1677       
1678    USE indices,                                                               &
1679        ONLY:  ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums,       &
1680               ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzb_diff_s_inner,         &
1681               nzb_s_inner, nzt, nzt_diff, rflags_invers
1682       
1683    USE kinds
1684   
1685    USE land_surface_model_mod,                                                &
1686        ONLY:   ghf_eb, land_surface, m_soil, nzb_soil, nzt_soil,              &
1687                qsws_eb, qsws_liq_eb, qsws_soil_eb, qsws_veg_eb, r_a, r_s,     &
1688                shf_eb, t_soil
1689
1690    USE netcdf_interface,                                                      &
1691        ONLY:  dots_rad, dots_soil
1692
1693    USE pegrid
1694   
1695    USE radiation_model_mod,                                                   &
1696        ONLY:  radiation, radiation_scheme, rad_net,                 &
1697               rad_lw_in, rad_lw_out, rad_sw_in, rad_sw_out
1698
1699#if defined ( __rrtmg )
1700    USE radiation_model_mod,                                                   &
1701        ONLY:  rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir, rad_lw_cs_hr,   &
1702               rad_lw_hr,  rad_sw_cs_hr, rad_sw_hr
1703#endif
1704
1705    USE statistics
1706
1707    IMPLICIT NONE
1708
1709    INTEGER(iwp) ::  i                   !<
1710    INTEGER(iwp) ::  j                   !<
1711    INTEGER(iwp) ::  k                   !<
1712    INTEGER(iwp) ::  k_surface_level     !<
1713    INTEGER(iwp) ::  nt                  !<
1714    INTEGER(iwp) ::  omp_get_thread_num  !<
1715    INTEGER(iwp) ::  sr                  !<
1716    INTEGER(iwp) ::  tn                  !<
1717   
1718    LOGICAL ::  first  !<
1719   
1720    REAL(wp) ::  dptdz_threshold  !<
1721    REAL(wp) ::  fac              !<
1722    REAL(wp) ::  height           !<
1723    REAL(wp) ::  pts              !<
1724    REAL(wp) ::  sums_l_eper      !<
1725    REAL(wp) ::  sums_l_etot      !<
1726    REAL(wp) ::  s1               !<
1727    REAL(wp) ::  s2               !<
1728    REAL(wp) ::  s3               !<
1729    REAL(wp) ::  s4               !<
1730    REAL(wp) ::  s5               !<
1731    REAL(wp) ::  s6               !<
1732    REAL(wp) ::  s7               !<
1733    REAL(wp) ::  ust              !<
1734    REAL(wp) ::  ust2             !<
1735    REAL(wp) ::  u2               !<
1736    REAL(wp) ::  vst              !<
1737    REAL(wp) ::  vst2             !<
1738    REAL(wp) ::  v2               !<
1739    REAL(wp) ::  w2               !<
1740    REAL(wp) ::  z_i(2)           !<
1741
1742    REAL(wp) ::  dptdz(nzb+1:nzt+1)    !<
1743    REAL(wp) ::  sums_ll(nzb:nzt+1,2)  !<
1744
1745    CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
1746
1747!
1748!-- To be on the safe side, check whether flow_statistics has already been
1749!-- called once after the current time step
1750    IF ( flow_statistics_called )  THEN
1751
1752       message_string = 'flow_statistics is called two times within one ' // &
1753                        'timestep'
1754       CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
1755
1756    ENDIF
1757
1758    !$acc data create( sums, sums_l )
1759    !$acc update device( hom )
1760
1761!
1762!-- Compute statistics for each (sub-)region
1763    DO  sr = 0, statistic_regions
1764
1765!
1766!--    Initialize (local) summation array
1767       sums_l = 0.0_wp
1768
1769!
1770!--    Store sums that have been computed in other subroutines in summation
1771!--    array
1772       sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
1773!--    WARNING: next line still has to be adjusted for OpenMP
1774       sums_l(:,21,0) = sums_wsts_bc_l(:,sr)  ! heat flux from advec_s_bc
1775       sums_l(nzb+9,pr_palm,0)  = sums_divold_l(sr)  ! old divergence from pres
1776       sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
1777
1778!
1779!--    When calcuating horizontally-averaged total (resolved- plus subgrid-
1780!--    scale) vertical fluxes and velocity variances by using commonly-
1781!--    applied Reynolds-based methods ( e.g. <w'pt'> = (w-<w>)*(pt-<pt>) )
1782!--    in combination with the 5th order advection scheme, pronounced
1783!--    artificial kinks could be observed in the vertical profiles near the
1784!--    surface. Please note: these kinks were not related to the model truth,
1785!--    i.e. these kinks are just related to an evaluation problem.   
1786!--    In order avoid these kinks, vertical fluxes and horizontal as well
1787!--    vertical velocity variances are calculated directly within the advection
1788!--    routines, according to the numerical discretization, to evaluate the
1789!--    statistical quantities as they will appear within the prognostic
1790!--    equations.
1791!--    Copy the turbulent quantities, evaluated in the advection routines to
1792!--    the local array sums_l() for further computations.
1793       IF ( ws_scheme_mom .AND. sr == 0 )  THEN
1794
1795!
1796!--       According to the Neumann bc for the horizontal velocity components,
1797!--       the corresponding fluxes has to satisfiy the same bc.
1798          IF ( ocean )  THEN
1799             sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
1800             sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)
1801          ENDIF
1802
1803          DO  i = 0, threads_per_task-1
1804!
1805!--          Swap the turbulent quantities evaluated in advec_ws.
1806             sums_l(:,13,i) = sums_wsus_ws_l(:,i)       ! w*u*
1807             sums_l(:,15,i) = sums_wsvs_ws_l(:,i)       ! w*v*
1808             sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2
1809             sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2
1810             sums_l(:,32,i) = sums_ws2_ws_l(:,i)        ! w*2
1811             sums_l(:,34,i) = sums_l(:,34,i) + 0.5_wp *                        &
1812                              ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +      &
1813                                sums_ws2_ws_l(:,i) )    ! e*
1814             DO  k = nzb, nzt
1815                sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5_wp * ( &
1816                                                      sums_us2_ws_l(k,i) +     &
1817                                                      sums_vs2_ws_l(k,i) +     &
1818                                                      sums_ws2_ws_l(k,i)     )
1819             ENDDO
1820          ENDDO
1821
1822       ENDIF
1823
1824       IF ( ws_scheme_sca .AND. sr == 0 )  THEN
1825
1826          DO  i = 0, threads_per_task-1
1827             sums_l(:,17,i) = sums_wspts_ws_l(:,i)      ! w*pt* from advec_s_ws
1828             IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
1829             IF ( humidity .OR. passive_scalar ) sums_l(:,49,i) =              &
1830                                                   sums_wsqs_ws_l(:,i) !w*q*
1831          ENDDO
1832
1833       ENDIF
1834!
1835!--    Horizontally averaged profiles of horizontal velocities and temperature.
1836!--    They must have been computed before, because they are already required
1837!--    for other horizontal averages.
1838       tn = 0
1839
1840       !$OMP PARALLEL PRIVATE( i, j, k, tn )
1841!$     tn = omp_get_thread_num()
1842
1843       !$acc update device( sums_l )
1844
1845       !$OMP DO
1846       !$acc parallel loop gang present( pt, rflags_invers, rmask, sums_l, u, v ) create( s1, s2, s3 )
1847       DO  k = nzb, nzt+1
1848          s1 = 0
1849          s2 = 0
1850          s3 = 0
1851          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
1852          DO  i = nxl, nxr
1853             DO  j =  nys, nyn
1854!
1855!--             k+1 is used in rflags since rflags is set 0 at surface points
1856                s1 = s1 + u(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1857                s2 = s2 + v(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1858                s3 = s3 + pt(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1859             ENDDO
1860          ENDDO
1861          sums_l(k,1,tn) = s1
1862          sums_l(k,2,tn) = s2
1863          sums_l(k,4,tn) = s3
1864       ENDDO
1865       !$acc end parallel loop
1866
1867!
1868!--    Horizontally averaged profile of salinity
1869       IF ( ocean )  THEN
1870          !$OMP DO
1871          !$acc parallel loop gang present( rflags_invers, rmask, sums_l, sa ) create( s1 )
1872          DO  k = nzb, nzt+1
1873             s1 = 0
1874             !$acc loop vector collapse( 2 ) reduction( +: s1 )
1875             DO  i = nxl, nxr
1876                DO  j =  nys, nyn
1877                   s1 = s1 + sa(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1878                ENDDO
1879             ENDDO
1880             sums_l(k,23,tn) = s1
1881          ENDDO
1882          !$acc end parallel loop
1883       ENDIF
1884
1885!
1886!--    Horizontally averaged profiles of virtual potential temperature,
1887!--    total water content, specific humidity and liquid water potential
1888!--    temperature
1889       IF ( humidity )  THEN
1890
1891          !$OMP DO
1892          !$acc parallel loop gang present( q, rflags_invers, rmask, sums_l, vpt ) create( s1, s2 )
1893          DO  k = nzb, nzt+1
1894             s1 = 0
1895             s2 = 0
1896             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
1897             DO  i = nxl, nxr
1898                DO  j =  nys, nyn
1899                   s1 = s1 + q(k,j,i)   * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1900                   s2 = s2 + vpt(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1901                ENDDO
1902             ENDDO
1903             sums_l(k,41,tn) = s1
1904             sums_l(k,44,tn) = s2
1905          ENDDO
1906          !$acc end parallel loop
1907
1908          IF ( cloud_physics )  THEN
1909             !$OMP DO
1910             !$acc parallel loop gang present( pt, q, ql, rflags_invers, rmask, sums_l ) create( s1, s2 )
1911             DO  k = nzb, nzt+1
1912                s1 = 0
1913                s2 = 0
1914                !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
1915                DO  i = nxl, nxr
1916                   DO  j =  nys, nyn
1917                      s1 = s1 + ( q(k,j,i) - ql(k,j,i) ) * &
1918                                rmask(j,i,sr) * rflags_invers(j,i,k+1)
1919                      s2 = s2 + ( pt(k,j,i) + l_d_cp*pt_d_t(k) * ql(k,j,i) ) * &
1920                                rmask(j,i,sr) * rflags_invers(j,i,k+1)
1921                   ENDDO
1922                ENDDO
1923                sums_l(k,42,tn) = s1
1924                sums_l(k,43,tn) = s2
1925             ENDDO
1926             !$acc end parallel loop
1927          ENDIF
1928       ENDIF
1929
1930!
1931!--    Horizontally averaged profiles of passive scalar
1932       IF ( passive_scalar )  THEN
1933          !$OMP DO
1934          !$acc parallel loop gang present( q, rflags_invers, rmask, sums_l ) create( s1 )
1935          DO  k = nzb, nzt+1
1936             s1 = 0
1937             !$acc loop vector collapse( 2 ) reduction( +: s1 )
1938             DO  i = nxl, nxr
1939                DO  j =  nys, nyn
1940                   s1 = s1 + q(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
1941                ENDDO
1942             ENDDO
1943             sums_l(k,41,tn) = s1
1944          ENDDO
1945          !$acc end parallel loop
1946       ENDIF
1947       !$OMP END PARALLEL
1948
1949!
1950!--    Summation of thread sums
1951       IF ( threads_per_task > 1 )  THEN
1952          DO  i = 1, threads_per_task-1
1953             !$acc parallel present( sums_l )
1954             sums_l(:,1,0) = sums_l(:,1,0) + sums_l(:,1,i)
1955             sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i)
1956             sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i)
1957             !$acc end parallel
1958             IF ( ocean )  THEN
1959                !$acc parallel present( sums_l )
1960                sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)
1961                !$acc end parallel
1962             ENDIF
1963             IF ( humidity )  THEN
1964                !$acc parallel present( sums_l )
1965                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
1966                sums_l(:,44,0) = sums_l(:,44,0) + sums_l(:,44,i)
1967                !$acc end parallel
1968                IF ( cloud_physics )  THEN
1969                   !$acc parallel present( sums_l )
1970                   sums_l(:,42,0) = sums_l(:,42,0) + sums_l(:,42,i)
1971                   sums_l(:,43,0) = sums_l(:,43,0) + sums_l(:,43,i)
1972                   !$acc end parallel
1973                ENDIF
1974             ENDIF
1975             IF ( passive_scalar )  THEN
1976                !$acc parallel present( sums_l )
1977                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
1978                !$acc end parallel
1979             ENDIF
1980          ENDDO
1981       ENDIF
1982
1983#if defined( __parallel )
1984!
1985!--    Compute total sum from local sums
1986       !$acc update host( sums_l )
1987       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1988       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL,  &
1989                           MPI_SUM, comm2d, ierr )
1990       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1991       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL,  &
1992                           MPI_SUM, comm2d, ierr )
1993       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1994       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL,  &
1995                           MPI_SUM, comm2d, ierr )
1996       IF ( ocean )  THEN
1997          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1998          CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb,       &
1999                              MPI_REAL, MPI_SUM, comm2d, ierr )
2000       ENDIF
2001       IF ( humidity ) THEN
2002          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2003          CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb,       &
2004                              MPI_REAL, MPI_SUM, comm2d, ierr )
2005          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2006          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb,       &
2007                              MPI_REAL, MPI_SUM, comm2d, ierr )
2008          IF ( cloud_physics ) THEN
2009             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2010             CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb,    &
2011                                 MPI_REAL, MPI_SUM, comm2d, ierr )
2012             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2013             CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb,    &
2014                                 MPI_REAL, MPI_SUM, comm2d, ierr )
2015          ENDIF
2016       ENDIF
2017
2018       IF ( passive_scalar )  THEN
2019          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
2020          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb,       &
2021                              MPI_REAL, MPI_SUM, comm2d, ierr )
2022       ENDIF
2023       !$acc update device( sums )
2024#else
2025       !$acc parallel present( sums, sums_l )
2026       sums(:,1) = sums_l(:,1,0)
2027       sums(:,2) = sums_l(:,2,0)
2028       sums(:,4) = sums_l(:,4,0)
2029       !$acc end parallel
2030       IF ( ocean )  THEN
2031          !$acc parallel present( sums, sums_l )
2032          sums(:,23) = sums_l(:,23,0)
2033          !$acc end parallel
2034       ENDIF
2035       IF ( humidity )  THEN
2036          !$acc parallel present( sums, sums_l )
2037          sums(:,44) = sums_l(:,44,0)
2038          sums(:,41) = sums_l(:,41,0)
2039          !$acc end parallel
2040          IF ( cloud_physics )  THEN
2041             !$acc parallel present( sums, sums_l )
2042             sums(:,42) = sums_l(:,42,0)
2043             sums(:,43) = sums_l(:,43,0)
2044             !$acc end parallel
2045          ENDIF
2046       ENDIF
2047       IF ( passive_scalar )  THEN
2048          !$acc parallel present( sums, sums_l )
2049          sums(:,41) = sums_l(:,41,0)
2050          !$acc end parallel
2051       ENDIF
2052#endif
2053
2054!
2055!--    Final values are obtained by division by the total number of grid points
2056!--    used for summation. After that store profiles.
2057       !$acc parallel present( hom, ngp_2dh, ngp_2dh_s_inner, sums )
2058       sums(:,1) = sums(:,1) / ngp_2dh(sr)
2059       sums(:,2) = sums(:,2) / ngp_2dh(sr)
2060       sums(:,4) = sums(:,4) / ngp_2dh_s_inner(:,sr)
2061       hom(:,1,1,sr) = sums(:,1)             ! u
2062       hom(:,1,2,sr) = sums(:,2)             ! v
2063       hom(:,1,4,sr) = sums(:,4)             ! pt
2064       !$acc end parallel
2065
2066!
2067!--    Salinity
2068       IF ( ocean )  THEN
2069          !$acc parallel present( hom, ngp_2dh_s_inner, sums )
2070          sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr)
2071          hom(:,1,23,sr) = sums(:,23)             ! sa
2072          !$acc end parallel
2073       ENDIF
2074
2075!
2076!--    Humidity and cloud parameters
2077       IF ( humidity ) THEN
2078          !$acc parallel present( hom, ngp_2dh_s_inner, sums )
2079          sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr)
2080          sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
2081          hom(:,1,44,sr) = sums(:,44)                ! vpt
2082          hom(:,1,41,sr) = sums(:,41)                ! qv (q)
2083          !$acc end parallel
2084          IF ( cloud_physics ) THEN
2085             !$acc parallel present( hom, ngp_2dh_s_inner, sums )
2086             sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr)
2087             sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr)
2088             hom(:,1,42,sr) = sums(:,42)             ! qv
2089             hom(:,1,43,sr) = sums(:,43)             ! pt
2090             !$acc end parallel
2091          ENDIF
2092       ENDIF
2093
2094!
2095!--    Passive scalar
2096       IF ( passive_scalar )  THEN
2097          !$acc parallel present( hom, ngp_2dh_s_inner, sums )
2098          sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
2099          hom(:,1,41,sr) = sums(:,41)                ! s (q)
2100          !$acc end parallel
2101       ENDIF
2102
2103!
2104!--    Horizontally averaged profiles of the remaining prognostic variables,
2105!--    variances, the total and the perturbation energy (single values in last
2106!--    column of sums_l) and some diagnostic quantities.
2107!--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
2108!--    ----  speaking the following k-loop would have to be split up and
2109!--          rearranged according to the staggered grid.
2110!--          However, this implies no error since staggered velocity components
2111!--          are zero at the walls and inside buildings.
2112       tn = 0
2113       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper,             &
2114       !$OMP                   sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2,  &
2115       !$OMP                   w2 )
2116!$     tn = omp_get_thread_num()
2117
2118       !$OMP DO
2119       !$acc parallel loop gang present( e, hom, kh, km, p, pt, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3, s4, s5, s6, s7 )
2120       DO  k = nzb, nzt+1
2121          s1 = 0
2122          s2 = 0
2123          s3 = 0
2124          s4 = 0
2125          s5 = 0
2126          s6 = 0
2127          s7 = 0
2128          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5, s6, s7 )
2129          DO  i = nxl, nxr
2130             DO  j =  nys, nyn
2131!
2132!--             Prognostic and diagnostic variables
2133                s1 = s1 + w(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2134                s2 = s2 + e(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2135                s3 = s3 + km(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2136                s4 = s4 + kh(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2137                s5 = s5 + p(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2138                s6 = s6 + ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr) * &
2139                          rflags_invers(j,i,k+1)
2140!
2141!--             Higher moments
2142!--             (Computation of the skewness of w further below)
2143                s7 = s7 + w(k,j,i)**3 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2144             ENDDO
2145          ENDDO
2146          sums_l(k,3,tn)  = s1
2147          sums_l(k,8,tn)  = s2
2148          sums_l(k,9,tn)  = s3
2149          sums_l(k,10,tn) = s4
2150          sums_l(k,40,tn) = s5
2151          sums_l(k,33,tn) = s6
2152          sums_l(k,38,tn) = s7
2153       ENDDO
2154       !$acc end parallel loop
2155
2156       IF ( humidity )  THEN
2157          !$OMP DO
2158          !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l ) create( s1 )
2159          DO  k = nzb, nzt+1
2160             s1 = 0
2161             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2162             DO  i = nxl, nxr
2163                DO  j =  nys, nyn
2164                   s1 = s1 + ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr) * &
2165                             rflags_invers(j,i,k+1)
2166                ENDDO
2167             ENDDO
2168             sums_l(k,70,tn) = s1
2169          ENDDO
2170          !$acc end parallel loop
2171       ENDIF
2172
2173!
2174!--    Total and perturbation energy for the total domain (being
2175!--    collected in the last column of sums_l).
2176       s1 = 0
2177       !$OMP DO
2178       !$acc parallel loop collapse(3) present( rflags_invers, rmask, u, v, w ) reduction(+:s1)
2179       DO  i = nxl, nxr
2180          DO  j =  nys, nyn
2181             DO  k = nzb, nzt+1
2182                s1 = s1 + 0.5_wp *                                             & 
2183                          ( u(k,j,i)**2 + v(k,j,i)**2 + w(k,j,i)**2 ) *        &
2184                          rmask(j,i,sr) * rflags_invers(j,i,k+1)
2185             ENDDO
2186          ENDDO
2187       ENDDO
2188       !$acc end parallel loop
2189       !$acc parallel present( sums_l )
2190       sums_l(nzb+4,pr_palm,tn) = s1
2191       !$acc end parallel
2192
2193       !$OMP DO
2194       !$acc parallel present( rmask, sums_l, us, usws, vsws, ts ) create( s1, s2, s3, s4 )
2195       s1 = 0
2196       s2 = 0
2197       s3 = 0
2198       s4 = 0
2199       !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4 )
2200       DO  i = nxl, nxr
2201          DO  j =  nys, nyn
2202!
2203!--          2D-arrays (being collected in the last column of sums_l)
2204             s1 = s1 + us(j,i)   * rmask(j,i,sr)
2205             s2 = s2 + usws(j,i) * rmask(j,i,sr)
2206             s3 = s3 + vsws(j,i) * rmask(j,i,sr)
2207             s4 = s4 + ts(j,i)   * rmask(j,i,sr)
2208          ENDDO
2209       ENDDO
2210       sums_l(nzb,pr_palm,tn)   = s1
2211       sums_l(nzb+1,pr_palm,tn) = s2
2212       sums_l(nzb+2,pr_palm,tn) = s3
2213       sums_l(nzb+3,pr_palm,tn) = s4
2214       !$acc end parallel
2215
2216       IF ( humidity )  THEN
2217          !$acc parallel present( qs, rmask, sums_l ) create( s1 )
2218          s1 = 0
2219          !$acc loop vector collapse( 2 ) reduction( +: s1 )
2220          DO  i = nxl, nxr
2221             DO  j =  nys, nyn
2222                s1 = s1 + qs(j,i) * rmask(j,i,sr)
2223             ENDDO
2224          ENDDO
2225          sums_l(nzb+12,pr_palm,tn) = s1
2226          !$acc end parallel
2227       ENDIF
2228
2229!
2230!--    Computation of statistics when ws-scheme is not used. Else these
2231!--    quantities are evaluated in the advection routines.
2232       IF ( .NOT. ws_scheme_mom  .OR.  sr /= 0 )  THEN
2233
2234          !$OMP DO
2235          !$acc parallel loop gang present( u, v, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3, s4, ust2, vst2, w2 )
2236          DO  k = nzb, nzt+1
2237             s1 = 0
2238             s2 = 0
2239             s3 = 0
2240             s4 = 0
2241             !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4 )
2242             DO  i = nxl, nxr
2243                DO  j =  nys, nyn
2244                   ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
2245                   vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
2246                   w2   = w(k,j,i)**2
2247
2248                   s1 = s1 + ust2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2249                   s2 = s2 + vst2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2250                   s3 = s3 + w2   * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2251!
2252!--                Perturbation energy
2253                   s4 = s4 + 0.5_wp * ( ust2 + vst2 + w2 ) * rmask(j,i,sr) *   &
2254                             rflags_invers(j,i,k+1)
2255                ENDDO
2256             ENDDO
2257             sums_l(k,30,tn) = s1
2258             sums_l(k,31,tn) = s2
2259             sums_l(k,32,tn) = s3
2260             sums_l(k,34,tn) = s4
2261          ENDDO
2262          !$acc end parallel loop
2263!
2264!--       Total perturbation TKE
2265          !$OMP DO
2266          !$acc parallel present( sums_l ) create( s1 )
2267          s1 = 0
2268          !$acc loop reduction( +: s1 )
2269          DO  k = nzb, nzt+1
2270             s1 = s1 + sums_l(k,34,tn)
2271          ENDDO
2272          sums_l(nzb+5,pr_palm,tn) = s1
2273          !$acc end parallel
2274
2275       ENDIF
2276
2277!
2278!--    Horizontally averaged profiles of the vertical fluxes
2279
2280!
2281!--    Subgridscale fluxes.
2282!--    WARNING: If a Prandtl-layer is used (k=nzb for flat terrain), the fluxes
2283!--    -------  should be calculated there in a different way. This is done
2284!--             in the next loop further below, where results from this loop are
2285!--             overwritten. However, THIS WORKS IN CASE OF FLAT TERRAIN ONLY!
2286!--             The non-flat case still has to be handled.
2287!--    NOTE: for simplicity, nzb_s_inner is used below, although
2288!--    ----  strictly speaking the following k-loop would have to be
2289!--          split up according to the staggered grid.
2290!--          However, this implies no error since staggered velocity
2291!--          components are zero at the walls and inside buildings.
2292       !$OMP DO
2293       !$acc parallel loop gang present( ddzu, kh, km, pt, u, v, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3 )
2294       DO  k = nzb, nzt_diff
2295          s1 = 0
2296          s2 = 0
2297          s3 = 0
2298          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
2299          DO  i = nxl, nxr
2300             DO  j = nys, nyn
2301
2302!
2303!--             Momentum flux w"u"
2304                s1 = s1 - 0.25_wp * (                                          &
2305                               km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
2306                                                           ) * (               &
2307                                   ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
2308                                 + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
2309                                                               )               &
2310                               * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2311!
2312!--             Momentum flux w"v"
2313                s2 = s2 - 0.25_wp * (                                          &
2314                               km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
2315                                                           ) * (               &
2316                                   ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
2317                                 + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
2318                                                               )               &
2319                               * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2320!
2321!--             Heat flux w"pt"
2322                s3 = s3 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )                 &
2323                                 * ( pt(k+1,j,i) - pt(k,j,i) )                 &
2324                                 * ddzu(k+1) * rmask(j,i,sr)                   &
2325                                 * rflags_invers(j,i,k+1)
2326             ENDDO
2327          ENDDO
2328          sums_l(k,12,tn) = s1
2329          sums_l(k,14,tn) = s2
2330          sums_l(k,16,tn) = s3
2331       ENDDO
2332       !$acc end parallel loop
2333
2334!
2335!--    Salinity flux w"sa"
2336       IF ( ocean )  THEN
2337          !$acc parallel loop gang present( ddzu, kh, sa, rflags_invers, rmask, sums_l ) create( s1 )
2338          DO  k = nzb, nzt_diff
2339             s1 = 0
2340             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2341             DO  i = nxl, nxr
2342                DO  j = nys, nyn
2343                   s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
2344                                    * ( sa(k+1,j,i) - sa(k,j,i) )              &
2345                                    * ddzu(k+1) * rmask(j,i,sr)                & 
2346                                    * rflags_invers(j,i,k+1)
2347                ENDDO
2348             ENDDO
2349             sums_l(k,65,tn) = s1
2350          ENDDO
2351          !$acc end parallel loop
2352       ENDIF
2353
2354!
2355!--    Buoyancy flux, water flux (humidity flux) w"q"
2356       IF ( humidity ) THEN
2357
2358          !$acc parallel loop gang present( ddzu, kh, q, vpt, rflags_invers, rmask, sums_l ) create( s1, s2 )
2359          DO  k = nzb, nzt_diff
2360             s1 = 0
2361             s2 = 0
2362             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2363             DO  i = nxl, nxr
2364                DO  j = nys, nyn
2365                   s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
2366                                    * ( vpt(k+1,j,i) - vpt(k,j,i) )            &
2367                                    * ddzu(k+1) * rmask(j,i,sr)                &
2368                                    * rflags_invers(j,i,k+1)
2369                   s2 = s2 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
2370                                    * ( q(k+1,j,i) - q(k,j,i) )                &
2371                                    * ddzu(k+1) * rmask(j,i,sr)                &
2372                                    * rflags_invers(j,i,k+1)
2373                ENDDO
2374             ENDDO
2375             sums_l(k,45,tn) = s1
2376             sums_l(k,48,tn) = s2
2377          ENDDO
2378          !$acc end parallel loop
2379
2380          IF ( cloud_physics ) THEN
2381
2382             !$acc parallel loop gang present( ddzu, kh, q, ql, rflags_invers, rmask, sums_l ) create( s1 )
2383             DO  k = nzb, nzt_diff
2384                s1 = 0
2385                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2386                DO  i = nxl, nxr
2387                   DO  j = nys, nyn
2388                      s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )           &
2389                                       *  ( ( q(k+1,j,i) - ql(k+1,j,i) )       &
2390                                          - ( q(k,j,i) - ql(k,j,i) ) )         &
2391                                       * ddzu(k+1) * rmask(j,i,sr)             & 
2392                                       * rflags_invers(j,i,k+1)
2393                   ENDDO
2394                ENDDO
2395                sums_l(k,51,tn) = s1
2396             ENDDO
2397             !$acc end parallel loop
2398
2399          ENDIF
2400
2401       ENDIF
2402!
2403!--    Passive scalar flux
2404       IF ( passive_scalar )  THEN
2405
2406          !$acc parallel loop gang present( ddzu, kh, q, rflags_invers, rmask, sums_l ) create( s1 )
2407          DO  k = nzb, nzt_diff
2408             s1 = 0
2409             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2410             DO  i = nxl, nxr
2411                DO  j = nys, nyn
2412                   s1 = s1 - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )              &
2413                                    * ( q(k+1,j,i) - q(k,j,i) )                &
2414                                    * ddzu(k+1) * rmask(j,i,sr)                &
2415                                    * rflags_invers(j,i,k+1)
2416                ENDDO
2417             ENDDO
2418             sums_l(k,48,tn) = s1
2419          ENDDO
2420          !$acc end parallel loop
2421
2422       ENDIF
2423
2424       IF ( use_surface_fluxes )  THEN
2425
2426          !$OMP DO
2427          !$acc parallel present( rmask, shf, sums_l, usws, vsws ) create( s1, s2, s3, s4, s5 )
2428          s1 = 0
2429          s2 = 0
2430          s3 = 0
2431          s4 = 0
2432          s5 = 0
2433          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5 )
2434          DO  i = nxl, nxr
2435             DO  j =  nys, nyn
2436!
2437!--             Subgridscale fluxes in the Prandtl layer
2438                s1 = s1 + usws(j,i) * rmask(j,i,sr)     ! w"u"
2439                s2 = s2 + vsws(j,i) * rmask(j,i,sr)     ! w"v"
2440                s3 = s3 + shf(j,i)  * rmask(j,i,sr)     ! w"pt"
2441                s4 = s4 + 0.0_wp * rmask(j,i,sr)        ! u"pt"
2442                s5 = s5 + 0.0_wp * rmask(j,i,sr)        ! v"pt"
2443             ENDDO
2444          ENDDO
2445          sums_l(nzb,12,tn) = s1
2446          sums_l(nzb,14,tn) = s2
2447          sums_l(nzb,16,tn) = s3
2448          sums_l(nzb,58,tn) = s4
2449          sums_l(nzb,61,tn) = s5
2450          !$acc end parallel
2451
2452          IF ( ocean )  THEN
2453
2454             !$OMP DO
2455             !$acc parallel present( rmask, saswsb, sums_l ) create( s1 )
2456             s1 = 0
2457             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2458             DO  i = nxl, nxr
2459                DO  j =  nys, nyn
2460                   s1 = s1 + saswsb(j,i) * rmask(j,i,sr)  ! w"sa"
2461                ENDDO
2462             ENDDO
2463             sums_l(nzb,65,tn) = s1
2464             !$acc end parallel
2465
2466          ENDIF
2467
2468          IF ( humidity )  THEN
2469
2470             !$OMP DO
2471             !$acc parallel present( pt, q, qsws, rmask, shf, sums_l ) create( s1, s2 )
2472             s1 = 0
2473             s2 = 0
2474             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2475             DO  i = nxl, nxr
2476                DO  j =  nys, nyn
2477                   s1 = s1 + qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
2478                   s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) * shf(j,i)    &
2479                               + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )
2480                ENDDO
2481             ENDDO
2482             sums_l(nzb,48,tn) = s1
2483             sums_l(nzb,45,tn) = s2
2484             !$acc end parallel
2485
2486             IF ( cloud_droplets )  THEN
2487
2488                !$OMP DO
2489                !$acc parallel present( pt, q, ql, qsws, rmask, shf, sums_l ) create( s1 )
2490                s1 = 0
2491                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2492                DO  i = nxl, nxr
2493                   DO  j =  nys, nyn
2494                      s1 = s1 + ( ( 1.0_wp +                                   &
2495                                    0.61_wp * q(nzb,j,i) - ql(nzb,j,i) ) *     &
2496                                  shf(j,i) + 0.61_wp * pt(nzb,j,i) * qsws(j,i) )
2497                   ENDDO
2498                ENDDO
2499                sums_l(nzb,45,tn) = s1
2500                !$acc end parallel
2501
2502             ENDIF
2503
2504             IF ( cloud_physics )  THEN
2505
2506                !$OMP DO
2507                !$acc parallel present( qsws, rmask, sums_l ) create( s1 )
2508                s1 = 0
2509                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2510                DO  i = nxl, nxr
2511                   DO  j =  nys, nyn
2512!
2513!--                   Formula does not work if ql(nzb) /= 0.0
2514                      s1 = s1 + qsws(j,i) * rmask(j,i,sr)   ! w"q" (w"qv")
2515                   ENDDO
2516                ENDDO
2517                sums_l(nzb,51,tn) = s1
2518                !$acc end parallel
2519
2520             ENDIF
2521
2522          ENDIF
2523
2524          IF ( passive_scalar )  THEN
2525
2526             !$OMP DO
2527             !$acc parallel present( qsws, rmask, sums_l ) create( s1 )
2528             s1 = 0
2529             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2530             DO  i = nxl, nxr
2531                DO  j =  nys, nyn
2532                   s1 = s1 + qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
2533                ENDDO
2534             ENDDO
2535             sums_l(nzb,48,tn) = s1
2536             !$acc end parallel
2537
2538          ENDIF
2539
2540       ENDIF
2541
2542!
2543!--    Subgridscale fluxes at the top surface
2544       IF ( use_top_fluxes )  THEN
2545
2546          !$OMP DO
2547          !$acc parallel present( rmask, sums_l, tswst, uswst, vswst ) create( s1, s2, s3, s4, s5 )
2548          s1 = 0
2549          s2 = 0
2550          s3 = 0
2551          s4 = 0
2552          s5 = 0
2553          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5 )
2554          DO  i = nxl, nxr
2555             DO  j =  nys, nyn
2556                s1 = s1 + uswst(j,i) * rmask(j,i,sr)    ! w"u"
2557                s2 = s2 + vswst(j,i) * rmask(j,i,sr)    ! w"v"
2558                s3 = s3 + tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
2559                s4 = s4 + 0.0_wp * rmask(j,i,sr)        ! u"pt"
2560                s5 = s5 + 0.0_wp * rmask(j,i,sr)        ! v"pt"
2561             ENDDO
2562          ENDDO
2563          sums_l(nzt:nzt+1,12,tn) = s1
2564          sums_l(nzt:nzt+1,14,tn) = s2
2565          sums_l(nzt:nzt+1,16,tn) = s3
2566          sums_l(nzt:nzt+1,58,tn) = s4
2567          sums_l(nzt:nzt+1,61,tn) = s5
2568          !$acc end parallel
2569
2570          IF ( ocean )  THEN
2571
2572             !$OMP DO
2573             !$acc parallel present( rmask, saswst, sums_l ) create( s1 )
2574             s1 = 0
2575             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2576             DO  i = nxl, nxr
2577                DO  j =  nys, nyn
2578                   s1 = s1 + saswst(j,i) * rmask(j,i,sr)  ! w"sa"
2579                ENDDO
2580             ENDDO
2581             sums_l(nzt,65,tn) = s1
2582             !$acc end parallel
2583
2584          ENDIF
2585
2586          IF ( humidity )  THEN
2587
2588             !$OMP DO
2589             !$acc parallel present( pt, q, qswst, rmask, tswst, sums_l ) create( s1, s2 )
2590             s1 = 0
2591             s2 = 0
2592             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2593             DO  i = nxl, nxr
2594                DO  j =  nys, nyn
2595                   s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
2596                   s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) * tswst(j,i) +&
2597                                 0.61_wp * pt(nzt,j,i) * qswst(j,i) )
2598                ENDDO
2599             ENDDO
2600             sums_l(nzt,48,tn) = s1
2601             sums_l(nzt,45,tn) = s2
2602             !$acc end parallel
2603
2604             IF ( cloud_droplets )  THEN
2605
2606                !$OMP DO
2607                !$acc parallel present( pt, q, ql, qswst, rmask, tswst, sums_l ) create( s1 )
2608                s1 = 0
2609                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2610                DO  i = nxl, nxr
2611                   DO  j =  nys, nyn
2612                      s1 = s1 + ( ( 1.0_wp +                                   &
2613                                    0.61_wp * q(nzt,j,i) - ql(nzt,j,i) ) *     &
2614                                  tswst(j,i) +                                 &
2615                                  0.61_wp * pt(nzt,j,i) * qswst(j,i) )
2616                   ENDDO
2617                ENDDO
2618                sums_l(nzt,45,tn) = s1
2619                !$acc end parallel
2620
2621             ENDIF
2622
2623             IF ( cloud_physics )  THEN
2624
2625                !$OMP DO
2626                !$acc parallel present( qswst, rmask, sums_l ) create( s1 )
2627                s1 = 0
2628                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2629                DO  i = nxl, nxr
2630                   DO  j =  nys, nyn
2631!
2632!--                   Formula does not work if ql(nzb) /= 0.0
2633                      s1 = s1 + qswst(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
2634                   ENDDO
2635                ENDDO
2636                sums_l(nzt,51,tn) = s1
2637                !$acc end parallel
2638
2639             ENDIF
2640
2641          ENDIF
2642
2643          IF ( passive_scalar )  THEN
2644
2645             !$OMP DO
2646             !$acc parallel present( qswst, rmask, sums_l ) create( s1 )
2647             s1 = 0
2648             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2649             DO  i = nxl, nxr
2650                DO  j =  nys, nyn
2651                   s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
2652                ENDDO
2653             ENDDO
2654             sums_l(nzt,48,tn) = s1
2655             !$acc end parallel
2656
2657          ENDIF
2658
2659       ENDIF
2660
2661!
2662!--    Resolved fluxes (can be computed for all horizontal points)
2663!--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
2664!--    ----  speaking the following k-loop would have to be split up and
2665!--          rearranged according to the staggered grid.
2666       !$acc parallel loop gang present( hom, pt, rflags_invers, rmask, sums_l, u, v, w ) create( s1, s2, s3 )
2667       DO  k = nzb, nzt_diff
2668          s1 = 0
2669          s2 = 0
2670          s3 = 0
2671          !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
2672          DO  i = nxl, nxr
2673             DO  j = nys, nyn
2674                ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) + &
2675                                 u(k+1,j,i) - hom(k+1,1,1,sr) )
2676                vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) + &
2677                                 v(k+1,j,i) - hom(k+1,1,2,sr) )
2678                pts = 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) + &
2679                                 pt(k+1,j,i) - hom(k+1,1,4,sr) )
2680!
2681!--             Higher moments
2682                s1 = s1 + pts * w(k,j,i)**2 * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2683                s2 = s2 + pts**2 * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2684!
2685!--             Energy flux w*e* (has to be adjusted?)
2686                s3 = s3 + w(k,j,i) * 0.5_wp * ( ust**2 + vst**2 + w(k,j,i)**2 )&
2687                                   * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2688             ENDDO
2689          ENDDO
2690          sums_l(k,35,tn) = s1
2691          sums_l(k,36,tn) = s2
2692          sums_l(k,37,tn) = s3
2693       ENDDO
2694       !$acc end parallel loop
2695
2696!
2697!--    Salinity flux and density (density does not belong to here,
2698!--    but so far there is no other suitable place to calculate)
2699       IF ( ocean )  THEN
2700
2701          IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
2702
2703             !$acc parallel loop gang present( hom, rflags_invers, rmask, sa, sums_l, w ) create( s1 )
2704             DO  k = nzb, nzt_diff
2705                s1 = 0
2706                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2707                DO  i = nxl, nxr
2708                   DO  j = nys, nyn
2709                      s1 = s1 + 0.5_wp * ( sa(k,j,i)   - hom(k,1,23,sr) +      &
2710                                           sa(k+1,j,i) - hom(k+1,1,23,sr) )    &
2711                                       * w(k,j,i) * rmask(j,i,sr)              & 
2712                                       * rflags_invers(j,i,k+1)
2713                   ENDDO
2714                ENDDO
2715                sums_l(k,66,tn) = s1
2716             ENDDO
2717             !$acc end parallel loop
2718
2719          ENDIF
2720
2721          !$acc parallel loop gang present( rflags_invers, rho, prho, rmask, sums_l ) create( s1, s2 )
2722          DO  k = nzb, nzt_diff
2723             s1 = 0
2724             s2 = 0
2725             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2726             DO  i = nxl, nxr
2727                DO  j = nys, nyn
2728                   s1 = s1 + rho(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2729                   s2 = s2 + prho(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2730                ENDDO
2731             ENDDO
2732             sums_l(k,64,tn) = s1
2733             sums_l(k,71,tn) = s2
2734          ENDDO
2735          !$acc end parallel loop
2736
2737       ENDIF
2738
2739!
2740!--    Buoyancy flux, water flux, humidity flux, liquid water
2741!--    content, rain drop concentration and rain water content
2742       IF ( humidity )  THEN
2743
2744          IF ( cloud_physics  .OR.  cloud_droplets )  THEN
2745
2746             !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, vpt, w ) create( s1 )
2747             DO  k = nzb, nzt_diff
2748                s1 = 0
2749                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2750                DO  i = nxl, nxr
2751                   DO  j = nys, nyn
2752                      s1 = s1 + 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +     &
2753                                           vpt(k+1,j,i) - hom(k+1,1,44,sr) ) * &
2754                                         w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2755                   ENDDO
2756                ENDDO
2757                sums_l(k,46,tn) = s1
2758             ENDDO
2759             !$acc end parallel loop
2760
2761             IF ( .NOT. cloud_droplets )  THEN
2762
2763                !$acc parallel loop gang present( hom, q, ql, rflags_invers, rmask, sums_l, w ) create( s1 )
2764                DO  k = nzb, nzt_diff
2765                   s1 = 0
2766                   !$acc loop vector collapse( 2 ) reduction( +: s1 )
2767                   DO  i = nxl, nxr
2768                      DO  j = nys, nyn
2769                         s1 = s1 + 0.5_wp * ( ( q(k,j,i)   - ql(k,j,i)   ) - hom(k,1,42,sr) +   &
2770                                              ( q(k+1,j,i) - ql(k+1,j,i) ) - hom(k+1,1,42,sr) ) &
2771                                          * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2772                      ENDDO
2773                   ENDDO
2774                   sums_l(k,52,tn) = s1
2775                ENDDO
2776                !$acc end parallel loop
2777
2778                IF ( microphysics_seifert )  THEN
2779
2780                   !$acc parallel loop gang present( qc, ql, rflags_invers, rmask, sums_l ) create( s1, s2 )
2781                   DO  k = nzb, nzt_diff
2782                      s1 = 0
2783                      !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2784                      DO  i = nxl, nxr
2785                         DO  j = nys, nyn
2786                            s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2787                            s2 = s2 + qc(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2788                         ENDDO
2789                      ENDDO
2790                      sums_l(k,54,tn) = s1
2791                      sums_l(k,75,tn) = s2
2792                   ENDDO
2793                   !$acc end parallel loop
2794
2795                   !$acc parallel loop gang present( nr, qr, prr, rflags_invers, rmask, sums_l ) create( s1, s2, s3 )
2796                   DO  k = nzb, nzt_diff
2797                      s1 = 0
2798                      s2 = 0
2799                      s3 = 0
2800                      !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 )
2801                      DO  i = nxl, nxr
2802                         DO  j = nys, nyn
2803                            s1 = s1 + nr(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2804                            s2 = s2 + qr(k,j,i)  * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2805                            s3 = s3 + prr(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2806                         ENDDO
2807                      ENDDO
2808                      sums_l(k,73,tn) = s1
2809                      sums_l(k,74,tn) = s2
2810                      sums_l(k,76,tn) = s3
2811                   ENDDO
2812                   !$acc end parallel loop
2813
2814                ELSE
2815
2816                   !$acc parallel loop gang present( ql, rflags_invers, rmask, sums_l ) create( s1 )
2817                   DO  k = nzb, nzt_diff
2818                      s1 = 0
2819                      !$acc loop vector collapse( 2 ) reduction( +: s1 )
2820                      DO  i = nxl, nxr
2821                         DO  j = nys, nyn
2822                            s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2823                         ENDDO
2824                      ENDDO
2825                      sums_l(k,54,tn) = s1
2826                   ENDDO
2827                   !$acc end parallel loop
2828
2829                ENDIF
2830
2831             ELSE
2832
2833                !$acc parallel loop gang present( ql, rflags_invers, rmask, sums_l ) create( s1 )
2834                DO  k = nzb, nzt_diff
2835                   s1 = 0
2836                   !$acc loop vector collapse( 2 ) reduction( +: s1 )
2837                   DO  i = nxl, nxr
2838                      DO  j = nys, nyn
2839                         s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2840                      ENDDO
2841                   ENDDO
2842                   sums_l(k,54,tn) = s1
2843                ENDDO
2844                !$acc end parallel loop
2845
2846             ENDIF
2847
2848          ELSE
2849
2850             IF( .NOT. ws_scheme_sca  .OR.  sr /= 0 )  THEN
2851
2852                !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, vpt, w ) create( s1 )
2853                DO  k = nzb, nzt_diff
2854                   s1 = 0
2855                   !$acc loop vector collapse( 2 ) reduction( +: s1 )
2856                   DO  i = nxl, nxr
2857                      DO  j = nys, nyn
2858                         s1 = s1 + 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +   &
2859                                              vpt(k+1,j,i) - hom(k+1,1,44,sr) ) &
2860                                          * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1)
2861                      ENDDO
2862                   ENDDO
2863                   sums_l(k,46,tn) = s1
2864                ENDDO
2865                !$acc end parallel loop
2866
2867             ELSEIF ( ws_scheme_sca  .AND.  sr == 0 )  THEN
2868
2869                !$acc parallel loop present( hom, sums_l )
2870                DO  k = nzb, nzt_diff
2871                   sums_l(k,46,tn) = ( 1.0_wp + 0.61_wp * hom(k,1,41,sr) ) * sums_l(k,17,tn) + &
2872                                                0.61_wp * hom(k,1,4,sr) * sums_l(k,49,tn)
2873                ENDDO
2874                !$acc end parallel loop
2875
2876             ENDIF
2877
2878          ENDIF
2879
2880       ENDIF
2881!
2882!--    Passive scalar flux
2883       IF ( passive_scalar  .AND.  ( .NOT. ws_scheme_sca  .OR.  sr /= 0 ) )  THEN
2884
2885          !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l, w ) create( s1 )
2886          DO  k = nzb, nzt_diff
2887             s1 = 0
2888             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2889             DO  i = nxl, nxr
2890                DO  j = nys, nyn
2891                   s1 = s1 + 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +          &
2892                                        q(k+1,j,i) - hom(k+1,1,41,sr) )        &
2893                                    * w(k,j,i) * rmask(j,i,sr)                 &
2894                                    * rflags_invers(j,i,k+1)
2895                ENDDO
2896             ENDDO
2897             sums_l(k,49,tn) = s1
2898          ENDDO
2899          !$acc end parallel loop
2900
2901       ENDIF
2902
2903!
2904!--    For speed optimization fluxes which have been computed in part directly
2905!--    inside the WS advection routines are treated seperatly
2906!--    Momentum fluxes first:
2907       IF ( .NOT. ws_scheme_mom  .OR.  sr /= 0  )  THEN
2908
2909          !$OMP DO
2910          !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, u, v, w ) create( s1, s2 )
2911          DO  k = nzb, nzt_diff
2912             s1 = 0
2913             s2 = 0
2914             !$acc loop vector collapse( 2 ) reduction( +: s1, s2 )
2915             DO  i = nxl, nxr
2916                DO  j = nys, nyn
2917                   ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +               &
2918                                    u(k+1,j,i) - hom(k+1,1,1,sr) )
2919                   vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +               &
2920                                    v(k+1,j,i) - hom(k+1,1,2,sr) )
2921!
2922!--                Momentum flux w*u*
2923                   s1 = s1 + 0.5_wp * ( w(k,j,i-1) + w(k,j,i) )                &
2924                                    * ust * rmask(j,i,sr)                      &
2925                                    * rflags_invers(j,i,k+1)
2926!
2927!--                Momentum flux w*v*
2928                   s2 = s2 + 0.5_wp * ( w(k,j-1,i) + w(k,j,i) )                &
2929                                    * vst * rmask(j,i,sr)                      & 
2930                                    * rflags_invers(j,i,k+1)
2931                ENDDO
2932             ENDDO
2933             sums_l(k,13,tn) = s1
2934             sums_l(k,15,tn) = s1
2935          ENDDO
2936          !$acc end parallel loop
2937
2938       ENDIF
2939
2940       IF ( .NOT. ws_scheme_sca  .OR.  sr /= 0 )  THEN
2941
2942          !$OMP DO
2943          !$acc parallel loop gang present( hom, pt, rflags_invers, rmask, sums_l, w ) create( s1 )
2944          DO  k = nzb, nzt_diff
2945             s1 = 0
2946             !$acc loop vector collapse( 2 ) reduction( +: s1 )
2947             DO  i = nxl, nxr
2948                DO  j = nys, nyn
2949!
2950!--                Vertical heat flux
2951                   s1 = s1 + 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) +          &
2952                                        pt(k+1,j,i) - hom(k+1,1,4,sr) )        &
2953                                    * w(k,j,i) * rmask(j,i,sr)                 & 
2954                                    * rflags_invers(j,i,k+1)
2955                ENDDO
2956             ENDDO
2957             sums_l(k,17,tn) = s1
2958          ENDDO
2959          !$acc end parallel loop
2960
2961          IF ( humidity )  THEN
2962
2963             !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l, w ) create( s1 )
2964             DO  k = nzb, nzt_diff
2965                s1 = 0
2966                !$acc loop vector collapse( 2 ) reduction( +: s1 )
2967                DO  i = nxl, nxr
2968                   DO  j = nys, nyn
2969                      s1 = s1 + 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +       &
2970                                           q(k+1,j,i) - hom(k+1,1,41,sr) )     &
2971                                       * w(k,j,i) * rmask(j,i,sr)              &
2972                                       * rflags_invers(j,i,k+1)
2973                   ENDDO
2974                ENDDO
2975                sums_l(k,49,tn) = s1
2976             ENDDO
2977             !$acc end parallel loop
2978
2979          ENDIF
2980
2981       ENDIF
2982
2983
2984!
2985!--    Density at top follows Neumann condition
2986       IF ( ocean )  THEN
2987          !$acc parallel present( sums_l )
2988          sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
2989          sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
2990          !$acc end parallel
2991       ENDIF
2992
2993!
2994!--    Divergence of vertical flux of resolved scale energy and pressure
2995!--    fluctuations as well as flux of pressure fluctuation itself (68).
2996!--    First calculate the products, then the divergence.
2997!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
2998       IF ( hom(nzb+1,2,55,0) /= 0.0_wp  .OR.  hom(nzb+1,2,68,0) /= 0.0_wp )  THEN
2999
3000          STOP '+++ openACC porting for vertical flux div of resolved scale TKE in flow_statistics is still missing'
3001          sums_ll = 0.0_wp  ! local array
3002
3003          !$OMP DO
3004          DO  i = nxl, nxr
3005             DO  j = nys, nyn
3006                DO  k = nzb_s_inner(j,i)+1, nzt
3007
3008                   sums_ll(k,1) = sums_ll(k,1) + 0.5_wp * w(k,j,i) * (         &
3009                  ( 0.25_wp * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1) )  &
3010                            - 0.5_wp * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) ) )**2&
3011                + ( 0.25_wp * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i) )  &
3012                            - 0.5_wp * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) ) )**2&
3013                + w(k,j,i)**2                                        )
3014
3015                   sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i)             &
3016                                               * ( p(k,j,i) + p(k+1,j,i) )
3017
3018                ENDDO
3019             ENDDO
3020          ENDDO
3021          sums_ll(0,1)     = 0.0_wp    ! because w is zero at the bottom
3022          sums_ll(nzt+1,1) = 0.0_wp
3023          sums_ll(0,2)     = 0.0_wp
3024          sums_ll(nzt+1,2) = 0.0_wp
3025
3026          DO  k = nzb+1, nzt
3027             sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
3028             sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
3029             sums_l(k,68,tn) = sums_ll(k,2)
3030          ENDDO
3031          sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
3032          sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
3033          sums_l(nzb,68,tn) = 0.0_wp    ! because w* = 0 at nzb
3034
3035       ENDIF
3036
3037!
3038!--    Divergence of vertical flux of SGS TKE and the flux itself (69)
3039       IF ( hom(nzb+1,2,57,0) /= 0.0_wp  .OR.  hom(nzb+1,2,69,0) /= 0.0_wp )  THEN
3040
3041          STOP '+++ openACC porting for vertical flux div of SGS TKE in flow_statistics is still missing'
3042          !$OMP DO
3043          DO  i = nxl, nxr
3044             DO  j = nys, nyn
3045                DO  k = nzb_s_inner(j,i)+1, nzt
3046
3047                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5_wp * (              &
3048                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
3049                 - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
3050                                                                ) * ddzw(k)
3051
3052                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5_wp * (              &
3053                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
3054                                                                )
3055
3056                ENDDO
3057             ENDDO
3058          ENDDO
3059          sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
3060          sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
3061
3062       ENDIF
3063
3064!
3065!--    Horizontal heat fluxes (subgrid, resolved, total).
3066!--    Do it only, if profiles shall be plotted.
3067       IF ( hom(nzb+1,2,58,0) /= 0.0_wp ) THEN
3068
3069          STOP '+++ openACC porting for horizontal flux calculation in flow_statistics is still missing'
3070          !$OMP DO
3071          DO  i = nxl, nxr
3072             DO  j = nys, nyn
3073                DO  k = nzb_s_inner(j,i)+1, nzt
3074!
3075!--                Subgrid horizontal heat fluxes u"pt", v"pt"
3076                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5_wp *                &
3077                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
3078                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
3079                                                 * ddx * rmask(j,i,sr)
3080                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
3081                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
3082                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
3083                                                 * ddy * rmask(j,i,sr)
3084!
3085!--                Resolved horizontal heat fluxes u*pt*, v*pt*
3086                   sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
3087                                     ( u(k,j,i) - hom(k,1,1,sr) ) * 0.5_wp *   &
3088                                     ( pt(k,j,i-1) - hom(k,1,4,sr) +           &
3089                                       pt(k,j,i)   - hom(k,1,4,sr) )
3090                   pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
3091                                    pt(k,j,i)   - hom(k,1,4,sr) )
3092                   sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
3093                                     ( v(k,j,i) - hom(k,1,2,sr) ) * 0.5_wp *   &
3094                                     ( pt(k,j-1,i) - hom(k,1,4,sr) +           &
3095                                       pt(k,j,i)   - hom(k,1,4,sr) )
3096                ENDDO
3097             ENDDO
3098          ENDDO
3099!
3100!--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
3101          sums_l(nzb,58,tn) = 0.0_wp
3102          sums_l(nzb,59,tn) = 0.0_wp
3103          sums_l(nzb,60,tn) = 0.0_wp
3104          sums_l(nzb,61,tn) = 0.0_wp
3105          sums_l(nzb,62,tn) = 0.0_wp
3106          sums_l(nzb,63,tn) = 0.0_wp
3107
3108       ENDIF
3109
3110!
3111!--    Collect current large scale advection and subsidence tendencies for
3112!--    data output
3113       IF ( large_scale_forcing  .AND.  ( simulated_time > 0.0_wp ) )  THEN
3114!
3115!--       Interpolation in time of LSF_DATA
3116          nt = 1
3117          DO WHILE ( simulated_time - dt_3d > time_vert(nt) )
3118             nt = nt + 1
3119          ENDDO
3120          IF ( simulated_time - dt_3d /= time_vert(nt) )  THEN
3121            nt = nt - 1
3122          ENDIF
3123
3124          fac = ( simulated_time - dt_3d - time_vert(nt) )                     &
3125                / ( time_vert(nt+1)-time_vert(nt) )
3126
3127
3128          DO  k = nzb, nzt
3129             sums_ls_l(k,0) = td_lsa_lpt(k,nt)                                 &
3130                              + fac * ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )
3131             sums_ls_l(k,1) = td_lsa_q(k,nt)                                   &
3132                              + fac * ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) )
3133          ENDDO
3134
3135          sums_ls_l(nzt+1,0) = sums_ls_l(nzt,0)
3136          sums_ls_l(nzt+1,1) = sums_ls_l(nzt,1)
3137
3138          IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
3139
3140             DO  k = nzb, nzt
3141                sums_ls_l(k,2) = td_sub_lpt(k,nt) + fac *                      &
3142                                 ( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )
3143                sums_ls_l(k,3) = td_sub_q(k,nt) + fac *                        &
3144                                 ( td_sub_q(k,nt+1) - td_sub_q(k,nt) )
3145             ENDDO
3146
3147             sums_ls_l(nzt+1,2) = sums_ls_l(nzt,2)
3148             sums_ls_l(nzt+1,3) = sums_ls_l(nzt,3)
3149
3150          ENDIF
3151
3152       ENDIF
3153
3154
3155       IF ( land_surface )  THEN
3156          !$OMP DO
3157          DO  i = nxl, nxr
3158             DO  j =  nys, nyn
3159                DO  k = nzb_soil, nzt_soil
3160                   sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil(k,j,i)         &
3161                                      * rmask(j,i,sr)
3162                   sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil(k,j,i)         &
3163                                      * rmask(j,i,sr)
3164                ENDDO
3165             ENDDO
3166          ENDDO
3167       ENDIF
3168
3169
3170       IF ( radiation .AND. radiation_scheme == 'rrtmg' )  THEN
3171          !$OMP DO
3172          DO  i = nxl, nxr
3173             DO  j =  nys, nyn
3174                DO  k = nzb_s_inner(j,i)+1, nzt+1
3175                   sums_l(k,102,tn)  = sums_l(k,102,tn)  + rad_lw_in(k,j,i)    &
3176                                       * rmask(j,i,sr)
3177                   sums_l(k,103,tn)  = sums_l(k,103,tn)  + rad_lw_out(k,j,i)   &
3178                                       * rmask(j,i,sr)
3179                   sums_l(k,104,tn)  = sums_l(k,104,tn)  + rad_sw_in(k,j,i)    &
3180                                       * rmask(j,i,sr)
3181                   sums_l(k,105,tn)  = sums_l(k,105,tn)  + rad_sw_out(k,j,i)   &
3182                                       * rmask(j,i,sr)
3183#if defined ( __rrtmg )
3184                   sums_l(k,106,tn)  = sums_l(k,106,tn)  + rad_lw_cs_hr(k,j,i) &
3185                                       * rmask(j,i,sr)
3186                   sums_l(k,107,tn)  = sums_l(k,107,tn)  + rad_lw_hr(k,j,i)    &
3187                                       * rmask(j,i,sr)
3188                   sums_l(k,108,tn)  = sums_l(k,108,tn)  + rad_sw_cs_hr(k,j,i) &
3189                                       * rmask(j,i,sr)
3190                   sums_l(k,109,tn)  = sums_l(k,109,tn)  + rad_sw_hr(k,j,i)    &
3191                                       * rmask(j,i,sr)
3192#endif
3193                ENDDO
3194             ENDDO
3195          ENDDO
3196       ENDIF
3197
3198!
3199!--    Calculate the user-defined profiles
3200       CALL user_statistics( 'profiles', sr, tn )
3201       !$OMP END PARALLEL
3202
3203!
3204!--    Summation of thread sums
3205       IF ( threads_per_task > 1 )  THEN
3206          STOP '+++ openACC porting for threads_per_task > 1 in flow_statistics is still missing'
3207          DO  i = 1, threads_per_task-1
3208             sums_l(:,3,0)          = sums_l(:,3,0) + sums_l(:,3,i)
3209             sums_l(:,4:40,0)       = sums_l(:,4:40,0) + sums_l(:,4:40,i)
3210             sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
3211                                      sums_l(:,45:pr_palm,i)
3212             IF ( max_pr_user > 0 )  THEN
3213                sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
3214                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
3215                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
3216             ENDIF
3217          ENDDO
3218       ENDIF
3219
3220       !$acc update host( hom, sums, sums_l )
3221
3222#if defined( __parallel )
3223
3224!
3225!--    Compute total sum from local sums
3226       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
3227       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, &
3228                           MPI_SUM, comm2d, ierr )
3229       IF ( large_scale_forcing )  THEN
3230          CALL MPI_ALLREDUCE( sums_ls_l(nzb,2), sums(nzb,83), ngp_sums_ls,     &
3231                              MPI_REAL, MPI_SUM, comm2d, ierr )
3232       ENDIF
3233#else
3234       sums = sums_l(:,:,0)
3235       IF ( large_scale_forcing )  THEN
3236          sums(:,81:88) = sums_ls_l
3237       ENDIF
3238#endif
3239
3240!
3241!--    Final values are obtained by division by the total number of grid points
3242!--    used for summation. After that store profiles.
3243!--    Check, if statistical regions do contain at least one grid point at the
3244!--    respective k-level, otherwise division by zero will lead to undefined
3245!--    values, which may cause e.g. problems with NetCDF output
3246!--    Profiles:
3247       DO  k = nzb, nzt+1
3248          sums(k,3)             = sums(k,3)             / ngp_2dh(sr)
3249          sums(k,12:22)         = sums(k,12:22)         / ngp_2dh(sr)
3250          sums(k,30:32)         = sums(k,30:32)         / ngp_2dh(sr)
3251          sums(k,35:39)         = sums(k,35:39)         / ngp_2dh(sr)
3252          sums(k,45:53)         = sums(k,45:53)         / ngp_2dh(sr)
3253          sums(k,55:63)         = sums(k,55:63)         / ngp_2dh(sr)
3254          sums(k,81:88)         = sums(k,81:88)         / ngp_2dh(sr)
3255          sums(k,89:114)        = sums(k,89:114)        / ngp_2dh(sr)
3256          IF ( ngp_2dh_s_inner(k,sr) /= 0 )  THEN
3257             sums(k,8:11)          = sums(k,8:11)          / ngp_2dh_s_inner(k,sr)
3258             sums(k,23:29)         = sums(k,23:29)         / ngp_2dh_s_inner(k,sr)
3259             sums(k,33:34)         = sums(k,33:34)         / ngp_2dh_s_inner(k,sr)
3260             sums(k,40)            = sums(k,40)            / ngp_2dh_s_inner(k,sr)
3261             sums(k,54)            = sums(k,54)            / ngp_2dh_s_inner(k,sr)
3262             sums(k,64)            = sums(k,64)            / ngp_2dh_s_inner(k,sr)
3263             sums(k,70:80)         = sums(k,70:80)         / ngp_2dh_s_inner(k,sr)
3264             sums(k,115:pr_palm-2) = sums(k,115:pr_palm-2) / ngp_2dh_s_inner(k,sr)
3265          ENDIF
3266       ENDDO
3267
3268!--    u* and so on
3269!--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
3270!--    size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
3271!--    above the topography, they are being divided by ngp_2dh(sr)
3272       sums(nzb:nzb+3,pr_palm)    = sums(nzb:nzb+3,pr_palm)    / &
3273                                    ngp_2dh(sr)
3274       sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)       / &    ! qs
3275                                    ngp_2dh(sr)
3276!--    eges, e*
3277       sums(nzb+4:nzb+5,pr_palm)  = sums(nzb+4:nzb+5,pr_palm)  / &
3278                                    ngp_3d(sr)
3279!--    Old and new divergence
3280       sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
3281                                    ngp_3d_inner(sr)
3282
3283!--    User-defined profiles
3284       IF ( max_pr_user > 0 )  THEN
3285          DO  k = nzb, nzt+1
3286             IF ( ngp_2dh_s_inner(k,sr) /= 0 )  THEN
3287                sums(k,pr_palm+1:pr_palm+max_pr_user) = &
3288                                       sums(k,pr_palm+1:pr_palm+max_pr_user) / &
3289                                       ngp_2dh_s_inner(k,sr)
3290             ENDIF
3291          ENDDO
3292       ENDIF
3293
3294!
3295!--    Collect horizontal average in hom.
3296!--    Compute deduced averages (e.g. total heat flux)
3297       hom(:,1,3,sr)  = sums(:,3)      ! w
3298       hom(:,1,8,sr)  = sums(:,8)      ! e     profiles 5-7 are initial profiles
3299       hom(:,1,9,sr)  = sums(:,9)      ! km
3300       hom(:,1,10,sr) = sums(:,10)     ! kh
3301       hom(:,1,11,sr) = sums(:,11)     ! l
3302       hom(:,1,12,sr) = sums(:,12)     ! w"u"
3303       hom(:,1,13,sr) = sums(:,13)     ! w*u*
3304       hom(:,1,14,sr) = sums(:,14)     ! w"v"
3305       hom(:,1,15,sr) = sums(:,15)     ! w*v*
3306       hom(:,1,16,sr) = sums(:,16)     ! w"pt"
3307       hom(:,1,17,sr) = sums(:,17)     ! w*pt*
3308       hom(:,1,18,sr) = sums(:,16) + sums(:,17)    ! wpt
3309       hom(:,1,19,sr) = sums(:,12) + sums(:,13)    ! wu
3310       hom(:,1,20,sr) = sums(:,14) + sums(:,15)    ! wv
3311       hom(:,1,21,sr) = sums(:,21)     ! w*pt*BC
3312       hom(:,1,22,sr) = sums(:,16) + sums(:,21)    ! wptBC
3313                                       ! profile 24 is initial profile (sa)
3314                                       ! profiles 25-29 left empty for initial
3315                                       ! profiles
3316       hom(:,1,30,sr) = sums(:,30)     ! u*2
3317       hom(:,1,31,sr) = sums(:,31)     ! v*2
3318       hom(:,1,32,sr) = sums(:,32)     ! w*2
3319       hom(:,1,33,sr) = sums(:,33)     ! pt*2
3320       hom(:,1,34,sr) = sums(:,34)     ! e*
3321       hom(:,1,35,sr) = sums(:,35)     ! w*2pt*
3322       hom(:,1,36,sr) = sums(:,36)     ! w*pt*2
3323       hom(:,1,37,sr) = sums(:,37)     ! w*e*
3324       hom(:,1,38,sr) = sums(:,38)     ! w*3
3325       hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20_wp )**1.5_wp   ! Sw
3326       hom(:,1,40,sr) = sums(:,40)     ! p
3327       hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
3328       hom(:,1,46,sr) = sums(:,46)     ! w*vpt*
3329       hom(:,1,47,sr) = sums(:,45) + sums(:,46)    ! wvpt
3330       hom(:,1,48,sr) = sums(:,48)     ! w"q" (w"qv")
3331       hom(:,1,49,sr) = sums(:,49)     ! w*q* (w*qv*)
3332       hom(:,1,50,sr) = sums(:,48) + sums(:,49)    ! wq (wqv)
3333       hom(:,1,51,sr) = sums(:,51)     ! w"qv"
3334       hom(:,1,52,sr) = sums(:,52)     ! w*qv*
3335       hom(:,1,53,sr) = sums(:,52) + sums(:,51)    ! wq (wqv)
3336       hom(:,1,54,sr) = sums(:,54)     ! ql
3337       hom(:,1,55,sr) = sums(:,55)     ! w*u*u*/dz
3338       hom(:,1,56,sr) = sums(:,56)     ! w*p*/dz
3339       hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho )/dz
3340       hom(:,1,58,sr) = sums(:,58)     ! u"pt"
3341       hom(:,1,59,sr) = sums(:,59)     ! u*pt*
3342       hom(:,1,60,sr) = sums(:,58) + sums(:,59)    ! upt_t
3343       hom(:,1,61,sr) = sums(:,61)     ! v"pt"
3344       hom(:,1,62,sr) = sums(:,62)     ! v*pt*
3345       hom(:,1,63,sr) = sums(:,61) + sums(:,62)    ! vpt_t
3346       hom(:,1,64,sr) = sums(:,64)     ! rho
3347       hom(:,1,65,sr) = sums(:,65)     ! w"sa"
3348       hom(:,1,66,sr) = sums(:,66)     ! w*sa*
3349       hom(:,1,67,sr) = sums(:,65) + sums(:,66)    ! wsa
3350       hom(:,1,68,sr) = sums(:,68)     ! w*p*
3351       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho
3352       hom(:,1,70,sr) = sums(:,70)     ! q*2
3353       hom(:,1,71,sr) = sums(:,71)     ! prho
3354       hom(:,1,72,sr) = hyp * 1E-4_wp     ! hyp in dbar
3355       hom(:,1,73,sr) = sums(:,73)     ! nr
3356       hom(:,1,74,sr) = sums(:,74)     ! qr
3357       hom(:,1,75,sr) = sums(:,75)     ! qc
3358       hom(:,1,76,sr) = sums(:,76)     ! prr (precipitation rate)
3359                                       ! 77 is initial density profile
3360       hom(:,1,78,sr) = ug             ! ug
3361       hom(:,1,79,sr) = vg             ! vg
3362       hom(:,1,80,sr) = w_subs         ! w_subs
3363
3364       IF ( large_scale_forcing )  THEN
3365          hom(:,1,81,sr) = sums_ls_l(:,0)          ! td_lsa_lpt
3366          hom(:,1,82,sr) = sums_ls_l(:,1)          ! td_lsa_q
3367          IF ( use_subsidence_tendencies )  THEN
3368             hom(:,1,83,sr) = sums_ls_l(:,2)       ! td_sub_lpt
3369             hom(:,1,84,sr) = sums_ls_l(:,3)       ! td_sub_q
3370          ELSE
3371             hom(:,1,83,sr) = sums(:,83)           ! td_sub_lpt
3372             hom(:,1,84,sr) = sums(:,84)           ! td_sub_q
3373          ENDIF
3374          hom(:,1,85,sr) = sums(:,85)              ! td_nud_lpt
3375          hom(:,1,86,sr) = sums(:,86)              ! td_nud_q
3376          hom(:,1,87,sr) = sums(:,87)              ! td_nud_u
3377          hom(:,1,88,sr) = sums(:,88)              ! td_nud_v
3378       END IF
3379
3380       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
3381                                       ! u*, w'u', w'v', t* (in last profile)
3382
3383       IF ( max_pr_user > 0 )  THEN    ! user-defined profiles
3384          hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
3385                               sums(:,pr_palm+1:pr_palm+max_pr_user)
3386       ENDIF
3387
3388!
3389!--    Determine the boundary layer height using two different schemes.
3390!--    First scheme: Starting from the Earth's (Ocean's) surface, look for the
3391!--    first relative minimum (maximum) of the total heat flux.
3392!--    The corresponding height is assumed as the boundary layer height, if it
3393!--    is less than 1.5 times the height where the heat flux becomes negative
3394!--    (positive) for the first time.
3395       z_i(1) = 0.0_wp
3396       first = .TRUE.
3397
3398       IF ( ocean )  THEN
3399          DO  k = nzt, nzb+1, -1
3400             IF (  first  .AND.  hom(k,1,18,sr) < -1.0E-8_wp )  THEN
3401                first = .FALSE.
3402                height = zw(k)
3403             ENDIF
3404             IF ( hom(k,1,18,sr) < -1.0E-8_wp  .AND.                           &
3405                  hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
3406                IF ( zw(k) < 1.5_wp * height )  THEN
3407                   z_i(1) = zw(k)
3408                ELSE
3409                   z_i(1) = height
3410                ENDIF
3411                EXIT
3412             ENDIF
3413          ENDDO
3414       ELSE
3415          DO  k = nzb, nzt-1
3416             IF ( first  .AND.  hom(k,1,18,sr) < -1.0E-8_wp )  THEN
3417                first = .FALSE.
3418                height = zw(k)
3419             ENDIF
3420             IF ( hom(k,1,18,sr) < -1.0E-8_wp  .AND.                           &
3421                  hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
3422                IF ( zw(k) < 1.5_wp * height )  THEN
3423                   z_i(1) = zw(k)
3424                ELSE
3425                   z_i(1) = height
3426                ENDIF
3427                EXIT
3428             ENDIF
3429          ENDDO
3430       ENDIF
3431
3432!
3433!--    Second scheme: Gradient scheme from Sullivan et al. (1998), modified
3434!--    by Uhlenbrock(2006). The boundary layer height is the height with the
3435!--    maximal local temperature gradient: starting from the second (the last
3436!--    but one) vertical gridpoint, the local gradient must be at least
3437!--    0.2K/100m and greater than the next four gradients.
3438!--    WARNING: The threshold value of 0.2K/100m must be adjusted for the
3439!--             ocean case!
3440       z_i(2) = 0.0_wp
3441       DO  k = nzb+1, nzt+1
3442          dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
3443       ENDDO
3444       dptdz_threshold = 0.2_wp / 100.0_wp
3445
3446       IF ( ocean )  THEN
3447          DO  k = nzt+1, nzb+5, -1
3448             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
3449                  dptdz(k) > dptdz(k-1)  .AND.  dptdz(k) > dptdz(k-2)  .AND.  &
3450                  dptdz(k) > dptdz(k-3)  .AND.  dptdz(k) > dptdz(k-4) )  THEN
3451                z_i(2) = zw(k-1)
3452                EXIT
3453             ENDIF
3454          ENDDO
3455       ELSE
3456          DO  k = nzb+1, nzt-3
3457             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
3458                  dptdz(k) > dptdz(k+1)  .AND.  dptdz(k) > dptdz(k+2)  .AND.  &
3459                  dptdz(k) > dptdz(k+3)  .AND.  dptdz(k) > dptdz(k+4) )  THEN
3460                z_i(2) = zw(k-1)
3461                EXIT
3462             ENDIF
3463          ENDDO
3464       ENDIF
3465
3466       hom(nzb+6,1,pr_palm,sr) = z_i(1)
3467       hom(nzb+7,1,pr_palm,sr) = z_i(2)
3468
3469!
3470!--    Determine vertical index which is nearest to the mean surface level
3471!--    height of the respective statistic region
3472       DO  k = nzb, nzt
3473          IF ( zw(k) >= mean_surface_level_height(sr) )  THEN
3474             k_surface_level = k
3475             EXIT
3476          ENDIF
3477       ENDDO
3478
3479!
3480!--    Computation of both the characteristic vertical velocity and
3481!--    the characteristic convective boundary layer temperature.
3482!--    The inversion height entering into the equation is defined with respect
3483!--    to the mean surface level height of the respective statistic region.
3484!--    The horizontal average at surface level index + 1 is input for the
3485!--    average temperature.
3486       IF ( hom(nzb,1,18,sr) > 1.0E-8_wp  .AND.  z_i(1) /= 0.0_wp )  THEN
3487          hom(nzb+8,1,pr_palm,sr) = &
3488             ( g / hom(k_surface_level+1,1,4,sr) * hom(k_surface_level,1,18,sr)&
3489             * ABS( z_i(1) - mean_surface_level_height(sr) ) )**0.333333333_wp
3490       ELSE
3491          hom(nzb+8,1,pr_palm,sr)  = 0.0_wp
3492       ENDIF
3493
3494!
3495!--    Collect the time series quantities
3496       ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)     ! E
3497       ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)     ! E*
3498       ts_value(3,sr) = dt_3d
3499       ts_value(4,sr) = hom(nzb,1,pr_palm,sr)       ! u*
3500       ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr)     ! th*
3501       ts_value(6,sr) = u_max
3502       ts_value(7,sr) = v_max
3503       ts_value(8,sr) = w_max
3504       ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr)    ! new divergence
3505       ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr)    ! old Divergence
3506       ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr)    ! z_i(1)
3507       ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr)    ! z_i(2)
3508       ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr)    ! w*
3509       ts_value(14,sr) = hom(nzb,1,16,sr)           ! w'pt'   at k=0
3510       ts_value(15,sr) = hom(nzb+1,1,16,sr)         ! w'pt'   at k=1
3511       ts_value(16,sr) = hom(nzb+1,1,18,sr)         ! wpt     at k=1
3512       ts_value(17,sr) = hom(nzb,1,4,sr)            ! pt(0)
3513       ts_value(18,sr) = hom(nzb+1,1,4,sr)          ! pt(zp)
3514       ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr)    ! u'w'    at k=0
3515       ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)    ! v'w'    at k=0
3516       ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
3517
3518       IF ( .NOT. neutral )  THEN
3519          ts_value(22,sr) = hom(nzb,1,114,sr)          ! L
3520       ELSE
3521          ts_value(22,sr) = 1.0E10_wp
3522       ENDIF
3523
3524       ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
3525
3526!
3527!--    Collect land surface model timeseries
3528       IF ( land_surface )  THEN
3529          ts_value(dots_soil  ,sr) = hom(nzb,1,93,sr)           ! ghf_eb
3530          ts_value(dots_soil+1,sr) = hom(nzb,1,94,sr)           ! shf_eb
3531          ts_value(dots_soil+2,sr) = hom(nzb,1,95,sr)           ! qsws_eb
3532          ts_value(dots_soil+3,sr) = hom(nzb,1,96,sr)           ! qsws_liq_eb
3533          ts_value(dots_soil+4,sr) = hom(nzb,1,97,sr)           ! qsws_soil_eb
3534          ts_value(dots_soil+5,sr) = hom(nzb,1,98,sr)           ! qsws_veg_eb
3535          ts_value(dots_soil+6,sr) = hom(nzb,1,99,sr)           ! r_a
3536          ts_value(dots_soil+7,sr) = hom(nzb,1,100,sr)          ! r_s
3537       ENDIF
3538!
3539!--    Collect radiation model timeseries
3540       IF ( radiation )  THEN
3541          ts_value(dots_rad,sr)   = hom(nzb,1,101,sr)          ! rad_net
3542          ts_value(dots_rad+1,sr) = hom(nzb,1,102,sr)          ! rad_lw_in
3543          ts_value(dots_rad+2,sr) = hom(nzb,1,103,sr)          ! rad_lw_out
3544          ts_value(dots_rad+3,sr) = hom(nzb,1,104,sr)          ! rad_sw_in
3545          ts_value(dots_rad+4,sr) = hom(nzb,1,105,sr)          ! rad_sw_out
3546
3547#if defined ( __rrtmg )
3548          IF ( radiation_scheme == 'rrtmg' )  THEN
3549             ts_value(dots_rad+5,sr) = hom(nzb,1,106,sr)          ! rrtm_aldif
3550             ts_value(dots_rad+6,sr) = hom(nzb,1,107,sr)          ! rrtm_aldir
3551             ts_value(dots_rad+7,sr) = hom(nzb,1,108,sr)          ! rrtm_asdif
3552             ts_value(dots_rad+8,sr) = hom(nzb,1,109,sr)          ! rrtm_asdir
3553          ENDIF
3554#endif
3555
3556       ENDIF
3557
3558!
3559!--    Calculate additional statistics provided by the user interface
3560       CALL user_statistics( 'time_series', sr, 0 )
3561
3562    ENDDO    ! loop of the subregions
3563
3564    !$acc end data
3565
3566!
3567!-- If required, sum up horizontal averages for subsequent time averaging
3568    IF ( do_sum )  THEN
3569       IF ( average_count_pr == 0 )  hom_sum = 0.0_wp
3570       hom_sum = hom_sum + hom(:,1,:,:)
3571       average_count_pr = average_count_pr + 1
3572       do_sum = .FALSE.
3573    ENDIF
3574
3575!
3576!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
3577!-- This flag is reset after each time step in time_integration.
3578    flow_statistics_called = .TRUE.
3579
3580    CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
3581
3582
3583 END SUBROUTINE flow_statistics
3584#endif
Note: See TracBrowser for help on using the repository browser.