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

Last change on this file since 2118 was 2118, checked in by raasch, 4 years ago

all OpenACC directives and related parts removed from the code

  • Property svn:keywords set to Id
File size: 77.2 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
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22! OpenACC version of subroutine removed
23!
24! Former revisions:
25! -----------------
26! $Id: flow_statistics.f90 2118 2017-01-17 16:38:49Z raasch $
27!
28! 2073 2016-11-30 14:34:05Z raasch
29! openmp bugfix: large scale forcing calculations cannot be executed thread
30! parallel
31!
32! 2037 2016-10-26 11:15:40Z knoop
33! Anelastic approximation implemented
34!
35! 2031 2016-10-21 15:11:58Z knoop
36! renamed variable rho to rho_ocean
37!
38! 2026 2016-10-18 10:27:02Z suehring
39! Bugfix, enable output of s*2.
40! Change, calculation of domain-averaged perturbation energy.
41! Some formatting adjustments.
42!
43! 2000 2016-08-20 18:09:15Z knoop
44! Forced header and separation lines into 80 columns
45!
46! 1976 2016-07-27 13:28:04Z maronga
47! Removed some unneeded __rrtmg preprocessor directives
48!
49! 1960 2016-07-12 16:34:24Z suehring
50! Separate humidity and passive scalar
51!
52! 1918 2016-05-27 14:35:57Z raasch
53! in case of Wicker-Skamarock scheme, calculate disturbance kinetic energy here,
54! if flow_statistics is called before the first initial time step
55!
56! 1853 2016-04-11 09:00:35Z maronga
57! Adjusted for use with radiation_scheme = constant
58!
59! 1849 2016-04-08 11:33:18Z hoffmann
60! prr moved to arrays_3d
61!
62! 1822 2016-04-07 07:49:42Z hoffmann
63! Output of bulk microphysics simplified.
64!
65! 1815 2016-04-06 13:49:59Z raasch
66! cpp-directives for intel openmp bug removed
67!
68! 1783 2016-03-06 18:36:17Z raasch
69! +module netcdf_interface
70!
71! 1747 2016-02-08 12:25:53Z raasch
72! small bugfixes for accelerator version
73!
74! 1738 2015-12-18 13:56:05Z raasch
75! bugfixes for calculations in statistical regions which do not contain grid
76! points in the lowest vertical levels, mean surface level height considered
77! in the calculation of the characteristic vertical velocity,
78! old upstream parts removed
79!
80! 1709 2015-11-04 14:47:01Z maronga
81! Updated output of Obukhov length
82!
83! 1691 2015-10-26 16:17:44Z maronga
84! Revised calculation of Obukhov length. Added output of radiative heating >
85! rates for RRTMG.
86!
87! 1682 2015-10-07 23:56:08Z knoop
88! Code annotations made doxygen readable
89!
90! 1658 2015-09-18 10:52:53Z raasch
91! bugfix: temporary reduction variables in the openacc branch are now
92! initialized to zero
93!
94! 1654 2015-09-17 09:20:17Z raasch
95! FORTRAN bugfix of r1652
96!
97! 1652 2015-09-17 08:12:24Z raasch
98! bugfix in calculation of energy production by turbulent transport of TKE
99!
100! 1593 2015-05-16 13:58:02Z raasch
101! FORTRAN errors removed from openacc branch
102!
103! 1585 2015-04-30 07:05:52Z maronga
104! Added output of timeseries and profiles for RRTMG
105!
106! 1571 2015-03-12 16:12:49Z maronga
107! Bugfix: output of rad_net and rad_sw_in
108!
109! 1567 2015-03-10 17:57:55Z suehring
110! Reverse modifications made for monotonic limiter.
111!
112! 1557 2015-03-05 16:43:04Z suehring
113! Adjustments for monotonic limiter
114!
115! 1555 2015-03-04 17:44:27Z maronga
116! Added output of r_a and r_s.
117!
118! 1551 2015-03-03 14:18:16Z maronga
119! Added suppport for land surface model and radiation model output.
120!
121! 1498 2014-12-03 14:09:51Z suehring
122! Comments added
123!
124! 1482 2014-10-18 12:34:45Z raasch
125! missing ngp_sums_ls added in accelerator version
126!
127! 1450 2014-08-21 07:31:51Z heinze
128! bugfix: calculate fac only for simulated_time >= 0.0
129!
130! 1396 2014-05-06 13:37:41Z raasch
131! bugfix: "copyin" replaced by "update device" in openacc-branch
132!
133! 1386 2014-05-05 13:55:30Z boeske
134! bugfix: simulated time before the last timestep is needed for the correct
135! calculation of the profiles of large scale forcing tendencies
136!
137! 1382 2014-04-30 12:15:41Z boeske
138! Renamed variables which store large scale forcing tendencies
139! pt_lsa -> td_lsa_lpt, pt_subs -> td_sub_lpt,
140! q_lsa  -> td_lsa_q,   q_subs  -> td_sub_q,
141! added Neumann boundary conditions for profile data output of large scale
142! advection and subsidence terms at nzt+1
143!
144! 1374 2014-04-25 12:55:07Z raasch
145! bugfix: syntax errors removed from openacc-branch
146! missing variables added to ONLY-lists
147!
148! 1365 2014-04-22 15:03:56Z boeske
149! Output of large scale advection, large scale subsidence and nudging tendencies
150! +sums_ls_l, ngp_sums_ls, use_subsidence_tendencies
151!
152! 1353 2014-04-08 15:21:23Z heinze
153! REAL constants provided with KIND-attribute
154!
155! 1322 2014-03-20 16:38:49Z raasch
156! REAL constants defined as wp-kind
157!
158! 1320 2014-03-20 08:40:49Z raasch
159! ONLY-attribute added to USE-statements,
160! kind-parameters added to all INTEGER and REAL declaration statements,
161! kinds are defined in new module kinds,
162! revision history before 2012 removed,
163! comment fields (!:) to be used for variable explanations added to
164! all variable declaration statements
165!
166! 1299 2014-03-06 13:15:21Z heinze
167! Output of large scale vertical velocity w_subs
168!
169! 1257 2013-11-08 15:18:40Z raasch
170! openacc "end parallel" replaced by "end parallel loop"
171!
172! 1241 2013-10-30 11:36:58Z heinze
173! Output of ug and vg
174!
175! 1221 2013-09-10 08:59:13Z raasch
176! ported for openACC in separate #else branch
177!
178! 1179 2013-06-14 05:57:58Z raasch
179! comment for profile 77 added
180!
181! 1115 2013-03-26 18:16:16Z hoffmann
182! ql is calculated by calc_liquid_water_content
183!
184! 1111 2013-03-08 23:54:10Z raasch
185! openACC directive added
186!
187! 1053 2012-11-13 17:11:03Z hoffmann
188! additions for two-moment cloud physics scheme:
189! +nr, qr, qc, prr
190!
191! 1036 2012-10-22 13:43:42Z raasch
192! code put under GPL (PALM 3.9)
193!
194! 1007 2012-09-19 14:30:36Z franke
195! Calculation of buoyancy flux for humidity in case of WS-scheme is now using
196! turbulent fluxes of WS-scheme
197! Bugfix: Calculation of subgridscale buoyancy flux for humidity and cloud
198! droplets at nzb and nzt added
199!
200! 801 2012-01-10 17:30:36Z suehring
201! Calculation of turbulent fluxes in advec_ws is now thread-safe.
202!
203! Revision 1.1  1997/08/11 06:15:17  raasch
204! Initial revision
205!
206!
207! Description:
208! ------------
209!> Compute average profiles and further average flow quantities for the different
210!> user-defined (sub-)regions. The region indexed 0 is the total model domain.
211!>
212!> @note For simplicity, nzb_s_inner and nzb_diff_s_inner are being used as a
213!>       lower vertical index for k-loops for all variables, although strictly
214!>       speaking the k-loops would have to be split up according to the staggered
215!>       grid. However, this implies no error since staggered velocity components
216!>       are zero at the walls and inside buildings.
217!------------------------------------------------------------------------------!
218 SUBROUTINE flow_statistics
219 
220
221    USE arrays_3d,                                                             &
222        ONLY:  ddzu, ddzw, e, heatflux_output_conversion, hyp, km, kh,         &
223               momentumflux_output_conversion, nr, ol, p, prho, prr, pt, q,    &
224               qc, ql, qr, qs, qsws, qswst, rho_air, rho_air_zw, rho_ocean, s, &
225               sa, ss, ssws, sswst, saswsb, saswst, shf, td_lsa_lpt, td_lsa_q, &
226               td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us, usws,    &
227               uswst, vsws, v, vg, vpt, vswst, w, w_subs,                      &
228               waterflux_output_conversion, zw
229       
230    USE cloud_parameters,                                                      &
231        ONLY:   l_d_cp, pt_d_t
232       
233    USE control_parameters,                                                    &
234        ONLY:   average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
235                dt_3d, g, humidity, kappa, large_scale_forcing,                &
236                large_scale_subsidence, max_pr_user, message_string, neutral,  &
237                microphysics_seifert, ocean, passive_scalar, simulated_time,   &
238                use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
239                ws_scheme_mom, ws_scheme_sca
240       
241    USE cpulog,                                                                &
242        ONLY:   cpu_log, log_point
243       
244    USE grid_variables,                                                        &
245        ONLY:   ddx, ddy
246       
247    USE indices,                                                               &
248        ONLY:   ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums,      &
249                ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzb_diff_s_inner,        &
250                nzb_s_inner, nzt, nzt_diff
251       
252    USE kinds
253   
254    USE land_surface_model_mod,                                                &
255        ONLY:   ghf_eb, land_surface, m_soil, nzb_soil, nzt_soil,              &
256                qsws_eb, qsws_liq_eb, qsws_soil_eb, qsws_veg_eb, r_a, r_s,     &
257                shf_eb, t_soil
258
259    USE netcdf_interface,                                                      &
260        ONLY:  dots_rad, dots_soil
261
262    USE pegrid
263
264    USE radiation_model_mod,                                                   &
265        ONLY:  radiation, radiation_scheme, rad_net,                           &
266               rad_lw_in, rad_lw_out, rad_lw_cs_hr, rad_lw_hr,                 &
267               rad_sw_in, rad_sw_out, rad_sw_cs_hr, rad_sw_hr
268
269#if defined ( __rrtmg )
270    USE radiation_model_mod,                                                   &
271        ONLY:  rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
272#endif
273 
274    USE statistics
275
276
277    IMPLICIT NONE
278
279    INTEGER(iwp) ::  i                   !<
280    INTEGER(iwp) ::  j                   !<
281    INTEGER(iwp) ::  k                   !<
282    INTEGER(iwp) ::  k_surface_level     !<
283    INTEGER(iwp) ::  nt                  !<
284    INTEGER(iwp) ::  omp_get_thread_num  !<
285    INTEGER(iwp) ::  sr                  !<
286    INTEGER(iwp) ::  tn                  !<
287   
288    LOGICAL ::  first  !<
289   
290    REAL(wp) ::  dptdz_threshold  !<
291    REAL(wp) ::  fac              !<
292    REAL(wp) ::  height           !<
293    REAL(wp) ::  pts              !<
294    REAL(wp) ::  sums_l_eper      !<
295    REAL(wp) ::  sums_l_etot      !<
296    REAL(wp) ::  ust              !<
297    REAL(wp) ::  ust2             !<
298    REAL(wp) ::  u2               !<
299    REAL(wp) ::  vst              !<
300    REAL(wp) ::  vst2             !<
301    REAL(wp) ::  v2               !<
302    REAL(wp) ::  w2               !<
303    REAL(wp) ::  z_i(2)           !<
304   
305    REAL(wp) ::  dptdz(nzb+1:nzt+1)    !<
306    REAL(wp) ::  sums_ll(nzb:nzt+1,2)  !<
307
308    CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
309
310
311!
312!-- To be on the safe side, check whether flow_statistics has already been
313!-- called once after the current time step
314    IF ( flow_statistics_called )  THEN
315
316       message_string = 'flow_statistics is called two times within one ' // &
317                        'timestep'
318       CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
319
320    ENDIF
321
322!
323!-- Compute statistics for each (sub-)region
324    DO  sr = 0, statistic_regions
325
326!
327!--    Initialize (local) summation array
328       sums_l = 0.0_wp
329
330!
331!--    Store sums that have been computed in other subroutines in summation
332!--    array
333       sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
334!--    WARNING: next line still has to be adjusted for OpenMP
335       sums_l(:,21,0) = sums_wsts_bc_l(:,sr) *                                 &
336                        heatflux_output_conversion  ! heat flux from advec_s_bc
337       sums_l(nzb+9,pr_palm,0)  = sums_divold_l(sr)  ! old divergence from pres
338       sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
339
340!
341!--    When calcuating horizontally-averaged total (resolved- plus subgrid-
342!--    scale) vertical fluxes and velocity variances by using commonly-
343!--    applied Reynolds-based methods ( e.g. <w'pt'> = (w-<w>)*(pt-<pt>) )
344!--    in combination with the 5th order advection scheme, pronounced
345!--    artificial kinks could be observed in the vertical profiles near the
346!--    surface. Please note: these kinks were not related to the model truth,
347!--    i.e. these kinks are just related to an evaluation problem.   
348!--    In order avoid these kinks, vertical fluxes and horizontal as well
349!--    vertical velocity variances are calculated directly within the advection
350!--    routines, according to the numerical discretization, to evaluate the
351!--    statistical quantities as they will appear within the prognostic
352!--    equations.
353!--    Copy the turbulent quantities, evaluated in the advection routines to
354!--    the local array sums_l() for further computations.
355       IF ( ws_scheme_mom .AND. sr == 0 )  THEN
356
357!
358!--       According to the Neumann bc for the horizontal velocity components,
359!--       the corresponding fluxes has to satisfiy the same bc.
360          IF ( ocean )  THEN
361             sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
362             sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)
363          ENDIF
364
365          DO  i = 0, threads_per_task-1
366!
367!--          Swap the turbulent quantities evaluated in advec_ws.
368             sums_l(:,13,i) = sums_wsus_ws_l(:,i)                              &
369                              * momentumflux_output_conversion ! w*u*
370             sums_l(:,15,i) = sums_wsvs_ws_l(:,i)                              &
371                              * momentumflux_output_conversion ! w*v*
372             sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2
373             sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2
374             sums_l(:,32,i) = sums_ws2_ws_l(:,i)        ! w*2
375             sums_l(:,34,i) = sums_l(:,34,i) + 0.5_wp *                        & 
376                              ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +      &
377                                sums_ws2_ws_l(:,i) )    ! e*
378          ENDDO
379
380       ENDIF
381
382       IF ( ws_scheme_sca .AND. sr == 0 )  THEN
383
384          DO  i = 0, threads_per_task-1
385             sums_l(:,17,i)                        = sums_wspts_ws_l(:,i)      &
386                                           * heatflux_output_conversion  ! w*pt*
387             IF ( ocean          ) sums_l(:,66,i)  = sums_wssas_ws_l(:,i) ! w*sa*
388             IF ( humidity       ) sums_l(:,49,i)  = sums_wsqs_ws_l(:,i)       &
389                                           * waterflux_output_conversion  ! w*q*
390             IF ( passive_scalar ) sums_l(:,116,i) = sums_wsss_ws_l(:,i)  ! w*s*
391          ENDDO
392
393       ENDIF
394!
395!--    Horizontally averaged profiles of horizontal velocities and temperature.
396!--    They must have been computed before, because they are already required
397!--    for other horizontal averages.
398       tn = 0
399
400       !$OMP PARALLEL PRIVATE( i, j, k, tn )
401!$     tn = omp_get_thread_num()
402
403       !$OMP DO
404       DO  i = nxl, nxr
405          DO  j =  nys, nyn
406             DO  k = nzb_s_inner(j,i), nzt+1
407                sums_l(k,1,tn)  = sums_l(k,1,tn)  + u(k,j,i)  * rmask(j,i,sr)
408                sums_l(k,2,tn)  = sums_l(k,2,tn)  + v(k,j,i)  * rmask(j,i,sr)
409                sums_l(k,4,tn)  = sums_l(k,4,tn)  + pt(k,j,i) * rmask(j,i,sr)
410             ENDDO
411          ENDDO
412       ENDDO
413
414!
415!--    Horizontally averaged profile of salinity
416       IF ( ocean )  THEN
417          !$OMP DO
418          DO  i = nxl, nxr
419             DO  j =  nys, nyn
420                DO  k = nzb_s_inner(j,i), nzt+1
421                   sums_l(k,23,tn)  = sums_l(k,23,tn) + &
422                                      sa(k,j,i) * rmask(j,i,sr)
423                ENDDO
424             ENDDO
425          ENDDO
426       ENDIF
427
428!
429!--    Horizontally averaged profiles of virtual potential temperature,
430!--    total water content, specific humidity and liquid water potential
431!--    temperature
432       IF ( humidity )  THEN
433          !$OMP DO
434          DO  i = nxl, nxr
435             DO  j =  nys, nyn
436                DO  k = nzb_s_inner(j,i), nzt+1
437                   sums_l(k,44,tn)  = sums_l(k,44,tn) + &
438                                      vpt(k,j,i) * rmask(j,i,sr)
439                   sums_l(k,41,tn)  = sums_l(k,41,tn) + &
440                                      q(k,j,i) * rmask(j,i,sr)
441                ENDDO
442             ENDDO
443          ENDDO
444          IF ( cloud_physics )  THEN
445             !$OMP DO
446             DO  i = nxl, nxr
447                DO  j =  nys, nyn
448                   DO  k = nzb_s_inner(j,i), nzt+1
449                      sums_l(k,42,tn) = sums_l(k,42,tn) + &
450                                      ( q(k,j,i) - ql(k,j,i) ) * rmask(j,i,sr)
451                      sums_l(k,43,tn) = sums_l(k,43,tn) + ( &
452                                      pt(k,j,i) + l_d_cp*pt_d_t(k) * ql(k,j,i) &
453                                                          ) * rmask(j,i,sr)
454                   ENDDO
455                ENDDO
456             ENDDO
457          ENDIF
458       ENDIF
459
460!
461!--    Horizontally averaged profiles of passive scalar
462       IF ( passive_scalar )  THEN
463          !$OMP DO
464          DO  i = nxl, nxr
465             DO  j =  nys, nyn
466                DO  k = nzb_s_inner(j,i), nzt+1
467                   sums_l(k,117,tn)  = sums_l(k,117,tn) + s(k,j,i) * rmask(j,i,sr)
468                ENDDO
469             ENDDO
470          ENDDO
471       ENDIF
472       !$OMP END PARALLEL
473!
474!--    Summation of thread sums
475       IF ( threads_per_task > 1 )  THEN
476          DO  i = 1, threads_per_task-1
477             sums_l(:,1,0) = sums_l(:,1,0) + sums_l(:,1,i)
478             sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i)
479             sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i)
480             IF ( ocean )  THEN
481                sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)
482             ENDIF
483             IF ( humidity )  THEN
484                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
485                sums_l(:,44,0) = sums_l(:,44,0) + sums_l(:,44,i)
486                IF ( cloud_physics )  THEN
487                   sums_l(:,42,0) = sums_l(:,42,0) + sums_l(:,42,i)
488                   sums_l(:,43,0) = sums_l(:,43,0) + sums_l(:,43,i)
489                ENDIF
490             ENDIF
491             IF ( passive_scalar )  THEN
492                sums_l(:,117,0) = sums_l(:,117,0) + sums_l(:,117,i)
493             ENDIF
494          ENDDO
495       ENDIF
496
497#if defined( __parallel )
498!
499!--    Compute total sum from local sums
500       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
501       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL,  &
502                           MPI_SUM, comm2d, ierr )
503       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
504       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL,  &
505                           MPI_SUM, comm2d, ierr )
506       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
507       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL,  &
508                           MPI_SUM, comm2d, ierr )
509       IF ( ocean )  THEN
510          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
511          CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb,       &
512                              MPI_REAL, MPI_SUM, comm2d, ierr )
513       ENDIF
514       IF ( humidity ) THEN
515          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
516          CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb,       &
517                              MPI_REAL, MPI_SUM, comm2d, ierr )
518          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
519          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb,       &
520                              MPI_REAL, MPI_SUM, comm2d, ierr )
521          IF ( cloud_physics ) THEN
522             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
523             CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb,    &
524                                 MPI_REAL, MPI_SUM, comm2d, ierr )
525             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
526             CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb,    &
527                                 MPI_REAL, MPI_SUM, comm2d, ierr )
528          ENDIF
529       ENDIF
530
531       IF ( passive_scalar )  THEN
532          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
533          CALL MPI_ALLREDUCE( sums_l(nzb,117,0), sums(nzb,117), nzt+2-nzb,       &
534                              MPI_REAL, MPI_SUM, comm2d, ierr )
535       ENDIF
536#else
537       sums(:,1) = sums_l(:,1,0)
538       sums(:,2) = sums_l(:,2,0)
539       sums(:,4) = sums_l(:,4,0)
540       IF ( ocean )  sums(:,23) = sums_l(:,23,0)
541       IF ( humidity ) THEN
542          sums(:,44) = sums_l(:,44,0)
543          sums(:,41) = sums_l(:,41,0)
544          IF ( cloud_physics ) THEN
545             sums(:,42) = sums_l(:,42,0)
546             sums(:,43) = sums_l(:,43,0)
547          ENDIF
548       ENDIF
549       IF ( passive_scalar )  sums(:,117) = sums_l(:,117,0)
550#endif
551
552!
553!--    Final values are obtained by division by the total number of grid points
554!--    used for summation. After that store profiles.
555       sums(:,1) = sums(:,1) / ngp_2dh(sr)
556       sums(:,2) = sums(:,2) / ngp_2dh(sr)
557       sums(:,4) = sums(:,4) / ngp_2dh_s_inner(:,sr)
558       hom(:,1,1,sr) = sums(:,1)             ! u
559       hom(:,1,2,sr) = sums(:,2)             ! v
560       hom(:,1,4,sr) = sums(:,4)             ! pt
561
562
563!
564!--    Salinity
565       IF ( ocean )  THEN
566          sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr)
567          hom(:,1,23,sr) = sums(:,23)             ! sa
568       ENDIF
569
570!
571!--    Humidity and cloud parameters
572       IF ( humidity ) THEN
573          sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr)
574          sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
575          hom(:,1,44,sr) = sums(:,44)             ! vpt
576          hom(:,1,41,sr) = sums(:,41)             ! qv (q)
577          IF ( cloud_physics ) THEN
578             sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr)
579             sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr)
580             hom(:,1,42,sr) = sums(:,42)             ! qv
581             hom(:,1,43,sr) = sums(:,43)             ! pt
582          ENDIF
583       ENDIF
584
585!
586!--    Passive scalar
587       IF ( passive_scalar )  hom(:,1,117,sr) = sums(:,117) /                  &
588            ngp_2dh_s_inner(:,sr)                    ! s
589
590!
591!--    Horizontally averaged profiles of the remaining prognostic variables,
592!--    variances, the total and the perturbation energy (single values in last
593!--    column of sums_l) and some diagnostic quantities.
594!--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
595!--    ----  speaking the following k-loop would have to be split up and
596!--          rearranged according to the staggered grid.
597!--          However, this implies no error since staggered velocity components
598!--          are zero at the walls and inside buildings.
599       tn = 0
600       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper,             &
601       !$OMP                   sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2,  &
602       !$OMP                   w2 )
603!$     tn = omp_get_thread_num()
604
605       !$OMP DO
606       DO  i = nxl, nxr
607          DO  j =  nys, nyn
608             sums_l_etot = 0.0_wp
609             DO  k = nzb_s_inner(j,i), nzt+1
610!
611!--             Prognostic and diagnostic variables
612                sums_l(k,3,tn)  = sums_l(k,3,tn)  + w(k,j,i)  * rmask(j,i,sr)
613                sums_l(k,8,tn)  = sums_l(k,8,tn)  + e(k,j,i)  * rmask(j,i,sr)
614                sums_l(k,9,tn)  = sums_l(k,9,tn)  + km(k,j,i) * rmask(j,i,sr)
615                sums_l(k,10,tn) = sums_l(k,10,tn) + kh(k,j,i) * rmask(j,i,sr)
616                sums_l(k,40,tn) = sums_l(k,40,tn) + p(k,j,i)
617
618                sums_l(k,33,tn) = sums_l(k,33,tn) + &
619                                  ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr)
620
621                IF ( humidity )  THEN
622                   sums_l(k,70,tn) = sums_l(k,70,tn) + &
623                                  ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr)
624                ENDIF
625                IF ( passive_scalar )  THEN
626                   sums_l(k,118,tn) = sums_l(k,118,tn) + &
627                                  ( s(k,j,i)-hom(k,1,117,sr) )**2 * rmask(j,i,sr)
628                ENDIF
629!
630!--             Higher moments
631!--             (Computation of the skewness of w further below)
632                sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i)**3 * rmask(j,i,sr)
633
634                sums_l_etot  = sums_l_etot + &
635                                        0.5_wp * ( u(k,j,i)**2 + v(k,j,i)**2 + &
636                                        w(k,j,i)**2 ) * rmask(j,i,sr)
637             ENDDO
638!
639!--          Total and perturbation energy for the total domain (being
640!--          collected in the last column of sums_l). Summation of these
641!--          quantities is seperated from the previous loop in order to
642!--          allow vectorization of that loop.
643             sums_l(nzb+4,pr_palm,tn) = sums_l(nzb+4,pr_palm,tn) + sums_l_etot
644!
645!--          2D-arrays (being collected in the last column of sums_l)
646             sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +               &
647                                        us(j,i)   * rmask(j,i,sr)
648             sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +             &
649                                        usws(j,i) * rmask(j,i,sr)
650             sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +             &
651                                        vsws(j,i) * rmask(j,i,sr)
652             sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +             &
653                                        ts(j,i)   * rmask(j,i,sr)
654             IF ( humidity )  THEN
655                sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +        &
656                                            qs(j,i)   * rmask(j,i,sr)
657             ENDIF
658             IF ( passive_scalar )  THEN
659                sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) +        &
660                                            ss(j,i)   * rmask(j,i,sr)
661             ENDIF
662          ENDDO
663       ENDDO
664
665!
666!--    Computation of statistics when ws-scheme is not used. Else these
667!--    quantities are evaluated in the advection routines.
668       IF ( .NOT. ws_scheme_mom .OR. sr /= 0 .OR. simulated_time == 0.0_wp )   &
669       THEN
670          !$OMP DO
671          DO  i = nxl, nxr
672             DO  j =  nys, nyn
673                DO  k = nzb_s_inner(j,i), nzt+1
674                   u2   = u(k,j,i)**2
675                   v2   = v(k,j,i)**2
676                   w2   = w(k,j,i)**2
677                   ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
678                   vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
679
680                   sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr)
681                   sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr)
682                   sums_l(k,32,tn) = sums_l(k,32,tn) + w2   * rmask(j,i,sr)
683!
684!--                Perturbation energy
685
686                   sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5_wp *                &
687                                  ( ust2 + vst2 + w2 ) * rmask(j,i,sr)
688                ENDDO
689             ENDDO
690          ENDDO
691       ENDIF
692!
693!--    Computaion of domain-averaged perturbation energy. Please note,
694!--    to prevent that perturbation energy is larger (even if only slightly)
695!--    than the total kinetic energy, calculation is based on deviations from
696!--    the horizontal mean, instead of spatial descretization of the advection
697!--    term.
698       !$OMP DO
699       DO  i = nxl, nxr
700          DO  j =  nys, nyn
701             DO  k = nzb_s_inner(j,i), nzt+1
702                w2   = w(k,j,i)**2
703                ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
704                vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
705                w2   = w(k,j,i)**2
706
707                sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn)            &
708                                 + 0.5_wp * ( ust2 + vst2 + w2 ) * rmask(j,i,sr)
709             ENDDO
710          ENDDO
711       ENDDO
712
713!
714!--    Horizontally averaged profiles of the vertical fluxes
715
716       !$OMP DO
717       DO  i = nxl, nxr
718          DO  j = nys, nyn
719!
720!--          Subgridscale fluxes (without Prandtl layer from k=nzb,
721!--          oterwise from k=nzb+1)
722!--          NOTE: for simplicity, nzb_diff_s_inner is used below, although
723!--          ----  strictly speaking the following k-loop would have to be
724!--                split up according to the staggered grid.
725!--                However, this implies no error since staggered velocity
726!--                components are zero at the walls and inside buildings.
727
728             DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
729!
730!--             Momentum flux w"u"
731                sums_l(k,12,tn) = sums_l(k,12,tn) - 0.25_wp * (                &
732                               km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
733                                                           ) * (               &
734                                   ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
735                                 + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
736                                                           ) * rmask(j,i,sr)   &
737                                         * rho_air_zw(k)                       &
738                                         * momentumflux_output_conversion(k)
739!
740!--             Momentum flux w"v"
741                sums_l(k,14,tn) = sums_l(k,14,tn) - 0.25_wp * (                &
742                               km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
743                                                           ) * (               &
744                                   ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
745                                 + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
746                                                           ) * rmask(j,i,sr)   &
747                                         * rho_air_zw(k)                       &
748                                         * momentumflux_output_conversion(k)
749!
750!--             Heat flux w"pt"
751                sums_l(k,16,tn) = sums_l(k,16,tn)                              &
752                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
753                                               * ( pt(k+1,j,i) - pt(k,j,i) )   &
754                                               * rho_air_zw(k)                 &
755                                               * heatflux_output_conversion(k) &
756                                               * ddzu(k+1) * rmask(j,i,sr)
757
758
759!
760!--             Salinity flux w"sa"
761                IF ( ocean )  THEN
762                   sums_l(k,65,tn) = sums_l(k,65,tn)                           &
763                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
764                                               * ( sa(k+1,j,i) - sa(k,j,i) )   &
765                                               * ddzu(k+1) * rmask(j,i,sr)
766                ENDIF
767
768!
769!--             Buoyancy flux, water flux (humidity flux) w"q"
770                IF ( humidity ) THEN
771                   sums_l(k,45,tn) = sums_l(k,45,tn)                           &
772                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
773                                               * ( vpt(k+1,j,i) - vpt(k,j,i) ) &
774                                               * rho_air_zw(k)                 &
775                                               * heatflux_output_conversion(k) &
776                                               * ddzu(k+1) * rmask(j,i,sr)
777                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
778                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
779                                               * ( q(k+1,j,i) - q(k,j,i) )     &
780                                               * rho_air_zw(k)                 &
781                                               * waterflux_output_conversion(k)&
782                                               * ddzu(k+1) * rmask(j,i,sr)
783
784                   IF ( cloud_physics ) THEN
785                      sums_l(k,51,tn) = sums_l(k,51,tn)                        &
786                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
787                                               * ( ( q(k+1,j,i) - ql(k+1,j,i) )&
788                                                - ( q(k,j,i) - ql(k,j,i) ) )   &
789                                               * rho_air_zw(k)                 &
790                                               * waterflux_output_conversion(k)&
791                                               * ddzu(k+1) * rmask(j,i,sr) 
792                   ENDIF
793                ENDIF
794
795!
796!--             Passive scalar flux
797                IF ( passive_scalar )  THEN
798                   sums_l(k,119,tn) = sums_l(k,119,tn)                         &
799                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
800                                                  * ( s(k+1,j,i) - s(k,j,i) )  &
801                                                  * ddzu(k+1) * rmask(j,i,sr)
802                ENDIF
803
804             ENDDO
805
806!
807!--          Subgridscale fluxes in the Prandtl layer
808             IF ( use_surface_fluxes )  THEN
809                sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
810                                    momentumflux_output_conversion(nzb) * &
811                                    usws(j,i) * rmask(j,i,sr)     ! w"u"
812                sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
813                                    momentumflux_output_conversion(nzb) * &
814                                    vsws(j,i) * rmask(j,i,sr)     ! w"v"
815                sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
816                                    heatflux_output_conversion(nzb) * &
817                                    shf(j,i)  * rmask(j,i,sr)     ! w"pt"
818                sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
819                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
820                sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
821                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
822                IF ( ocean )  THEN
823                   sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
824                                       saswsb(j,i) * rmask(j,i,sr)  ! w"sa"
825                ENDIF
826                IF ( humidity )  THEN
827                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
828                                       waterflux_output_conversion(nzb) *      &
829                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
830                   sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                   &
831                                       ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) *     &
832                                       shf(j,i) + 0.61_wp * pt(nzb,j,i) *      &
833                                                  qsws(j,i) )                  &
834                                       * heatflux_output_conversion(nzb)
835                   IF ( cloud_droplets )  THEN
836                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                &
837                                         ( 1.0_wp + 0.61_wp * q(nzb,j,i) -     &
838                                           ql(nzb,j,i) ) * shf(j,i) +          &
839                                           0.61_wp * pt(nzb,j,i) * qsws(j,i) ) &
840                                          * heatflux_output_conversion(nzb)
841                   ENDIF
842                   IF ( cloud_physics )  THEN
843!
844!--                   Formula does not work if ql(nzb) /= 0.0
845                      sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                  &
846                                          waterflux_output_conversion(nzb) *   &
847                                          qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
848                   ENDIF
849                ENDIF
850                IF ( passive_scalar )  THEN
851                   sums_l(nzb,119,tn) = sums_l(nzb,119,tn) +                     &
852                                        ssws(j,i) * rmask(j,i,sr) ! w"s"
853                ENDIF
854             ENDIF
855
856             IF ( .NOT. neutral )  THEN
857                sums_l(nzb,114,tn) = sums_l(nzb,114,tn) +                      &
858                                    ol(j,i)  * rmask(j,i,sr) ! L
859             ENDIF
860
861
862             IF ( land_surface )  THEN
863                sums_l(nzb,93,tn)  = sums_l(nzb,93,tn) + ghf_eb(j,i)
864                sums_l(nzb,94,tn)  = sums_l(nzb,94,tn) + shf_eb(j,i)
865                sums_l(nzb,95,tn)  = sums_l(nzb,95,tn) + qsws_eb(j,i)
866                sums_l(nzb,96,tn)  = sums_l(nzb,96,tn) + qsws_liq_eb(j,i)
867                sums_l(nzb,97,tn)  = sums_l(nzb,97,tn) + qsws_soil_eb(j,i)
868                sums_l(nzb,98,tn)  = sums_l(nzb,98,tn) + qsws_veg_eb(j,i)
869                sums_l(nzb,99,tn)  = sums_l(nzb,99,tn) + r_a(j,i)
870                sums_l(nzb,100,tn) = sums_l(nzb,100,tn)+ r_s(j,i)
871             ENDIF
872
873             IF ( radiation  .AND.  radiation_scheme /= 'constant' )  THEN
874                sums_l(nzb,101,tn)  = sums_l(nzb,101,tn)  + rad_net(j,i)
875                sums_l(nzb,102,tn)  = sums_l(nzb,102,tn)  + rad_lw_in(nzb,j,i)
876                sums_l(nzb,103,tn)  = sums_l(nzb,103,tn)  + rad_lw_out(nzb,j,i)
877                sums_l(nzb,104,tn)  = sums_l(nzb,104,tn)  + rad_sw_in(nzb,j,i)
878                sums_l(nzb,105,tn)  = sums_l(nzb,105,tn)  + rad_sw_out(nzb,j,i)
879
880#if defined ( __rrtmg )
881                IF ( radiation_scheme == 'rrtmg' )  THEN
882                   sums_l(nzb,110,tn)  = sums_l(nzb,110,tn)  + rrtm_aldif(0,j,i)
883                   sums_l(nzb,111,tn)  = sums_l(nzb,111,tn)  + rrtm_aldir(0,j,i)
884                   sums_l(nzb,112,tn)  = sums_l(nzb,112,tn)  + rrtm_asdif(0,j,i)
885                   sums_l(nzb,113,tn)  = sums_l(nzb,113,tn)  + rrtm_asdir(0,j,i)
886                ENDIF
887#endif
888             ENDIF
889!
890!--          Subgridscale fluxes at the top surface
891             IF ( use_top_fluxes )  THEN
892                sums_l(nzt:nzt+1,12,tn) = sums_l(nzt:nzt+1,12,tn) + &
893                                    momentumflux_output_conversion(nzt:nzt+1) * &
894                                    uswst(j,i) * rmask(j,i,sr)    ! w"u"
895                sums_l(nzt:nzt+1,14,tn) = sums_l(nzt:nzt+1,14,tn) + &
896                                    momentumflux_output_conversion(nzt:nzt+1) * &
897                                    vswst(j,i) * rmask(j,i,sr)    ! w"v"
898                sums_l(nzt:nzt+1,16,tn) = sums_l(nzt:nzt+1,16,tn) + &
899                                    heatflux_output_conversion(nzt:nzt+1) * &
900                                    tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
901                sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + &
902                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
903                sums_l(nzt:nzt+1,61,tn) = sums_l(nzt:nzt+1,61,tn) + &
904                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
905
906                IF ( ocean )  THEN
907                   sums_l(nzt,65,tn) = sums_l(nzt,65,tn) + &
908                                       saswst(j,i) * rmask(j,i,sr)  ! w"sa"
909                ENDIF
910                IF ( humidity )  THEN
911                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) +                     &
912                                       waterflux_output_conversion(nzt) *      &
913                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
914                   sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                   &
915                                       ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) *     &
916                                       tswst(j,i) + 0.61_wp * pt(nzt,j,i) *    &
917                                                             qswst(j,i) )      &
918                                       * heatflux_output_conversion(nzt)
919                   IF ( cloud_droplets )  THEN
920                      sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                &
921                                          ( 1.0_wp + 0.61_wp * q(nzt,j,i) -    &
922                                            ql(nzt,j,i) ) * tswst(j,i) +       &
923                                           0.61_wp * pt(nzt,j,i) * qswst(j,i) )&
924                                           * heatflux_output_conversion(nzt)
925                   ENDIF
926                   IF ( cloud_physics )  THEN
927!
928!--                   Formula does not work if ql(nzb) /= 0.0
929                      sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + &   ! w"q" (w"qv")
930                                          waterflux_output_conversion(nzt) *   &
931                                          qswst(j,i) * rmask(j,i,sr)
932                   ENDIF
933                ENDIF
934                IF ( passive_scalar )  THEN
935                   sums_l(nzt,119,tn) = sums_l(nzt,119,tn) + &
936                                        sswst(j,i) * rmask(j,i,sr) ! w"s"
937                ENDIF
938             ENDIF
939
940!
941!--          Resolved fluxes (can be computed for all horizontal points)
942!--          NOTE: for simplicity, nzb_s_inner is used below, although strictly
943!--          ----  speaking the following k-loop would have to be split up and
944!--                rearranged according to the staggered grid.
945             DO  k = nzb_s_inner(j,i), nzt
946                ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +                  &
947                                 u(k+1,j,i) - hom(k+1,1,1,sr) )
948                vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +                  &
949                                 v(k+1,j,i) - hom(k+1,1,2,sr) )
950                pts = 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) +                 &
951                                 pt(k+1,j,i) - hom(k+1,1,4,sr) )
952
953!--             Higher moments
954                sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 *        &
955                                                    rmask(j,i,sr)
956                sums_l(k,36,tn) = sums_l(k,36,tn) + pts**2 * w(k,j,i) *        &
957                                                    rmask(j,i,sr)
958
959!
960!--             Salinity flux and density (density does not belong to here,
961!--             but so far there is no other suitable place to calculate)
962                IF ( ocean )  THEN
963                   IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
964                      pts = 0.5_wp * ( sa(k,j,i)   - hom(k,1,23,sr) +          &
965                                       sa(k+1,j,i) - hom(k+1,1,23,sr) )
966                      sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) *     &
967                                        rmask(j,i,sr)
968                   ENDIF
969                   sums_l(k,64,tn) = sums_l(k,64,tn) + rho_ocean(k,j,i) *            &
970                                                       rmask(j,i,sr)
971                   sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) *           &
972                                                       rmask(j,i,sr)
973                ENDIF
974
975!
976!--             Buoyancy flux, water flux, humidity flux, liquid water
977!--             content, rain drop concentration and rain water content
978                IF ( humidity )  THEN
979                   IF ( cloud_physics .OR. cloud_droplets )  THEN
980                      pts = 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +         &
981                                    vpt(k+1,j,i) - hom(k+1,1,44,sr) )
982                      sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *     &
983                                               heatflux_output_conversion(k) * &
984                                                          rmask(j,i,sr)
985                      sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * rmask(j,i,sr)
986
987                      IF ( .NOT. cloud_droplets )  THEN
988                         pts = 0.5_wp *                                        &
989                              ( ( q(k,j,i) - ql(k,j,i) ) -                     &
990                              hom(k,1,42,sr) +                                 &
991                              ( q(k+1,j,i) - ql(k+1,j,i) ) -                   &
992                              hom(k+1,1,42,sr) )
993                         sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) *  &
994                                             waterflux_output_conversion(k) *  &
995                                                             rmask(j,i,sr)
996                         sums_l(k,75,tn) = sums_l(k,75,tn) + qc(k,j,i) *       &
997                                                             rmask(j,i,sr)
998                         sums_l(k,76,tn) = sums_l(k,76,tn) + prr(k,j,i) *      &
999                                                             rmask(j,i,sr)
1000                         IF ( microphysics_seifert )  THEN
1001                            sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) *    &
1002                                                                rmask(j,i,sr)
1003                            sums_l(k,74,tn) = sums_l(k,74,tn) + qr(k,j,i) *    &
1004                                                                rmask(j,i,sr)
1005                         ENDIF
1006                      ENDIF
1007
1008                   ELSE
1009                      IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
1010                         pts = 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +      &
1011                                          vpt(k+1,j,i) - hom(k+1,1,44,sr) )
1012                         sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *  &
1013                                              heatflux_output_conversion(k) *  &
1014                                                             rmask(j,i,sr)
1015                      ELSE IF ( ws_scheme_sca .AND. sr == 0 )  THEN
1016                         sums_l(k,46,tn) = ( ( 1.0_wp + 0.61_wp *              & 
1017                                               hom(k,1,41,sr) ) *              &
1018                                             sums_l(k,17,tn) +                 &
1019                                             0.61_wp * hom(k,1,4,sr) *         &
1020                                             sums_l(k,49,tn)                   &
1021                                           ) * heatflux_output_conversion(k)
1022                      END IF
1023                   END IF
1024                ENDIF
1025!
1026!--             Passive scalar flux
1027                IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca                &
1028                     .OR. sr /= 0 ) )  THEN
1029                   pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +             &
1030                                    s(k+1,j,i) - hom(k+1,1,117,sr) )
1031                   sums_l(k,116,tn) = sums_l(k,116,tn) + pts * w(k,j,i) *      &
1032                                                       rmask(j,i,sr)
1033                ENDIF
1034
1035!
1036!--             Energy flux w*e*
1037!--             has to be adjusted
1038                sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5_wp *        &
1039                                             ( ust**2 + vst**2 + w(k,j,i)**2 ) &
1040                                           * momentumflux_output_conversion(k) &
1041                                             * rmask(j,i,sr)
1042             ENDDO
1043          ENDDO
1044       ENDDO
1045!
1046!--    For speed optimization fluxes which have been computed in part directly
1047!--    inside the WS advection routines are treated seperatly
1048!--    Momentum fluxes first:
1049       IF ( .NOT. ws_scheme_mom .OR. sr /= 0  )  THEN
1050         !$OMP DO
1051         DO  i = nxl, nxr
1052            DO  j = nys, nyn
1053               DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
1054                  ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +                &
1055                                   u(k+1,j,i) - hom(k+1,1,1,sr) )
1056                  vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +                &
1057                                   v(k+1,j,i) - hom(k+1,1,2,sr) )
1058!
1059!--               Momentum flux w*u*
1060                  sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5_wp *                 &
1061                                                    ( w(k,j,i-1) + w(k,j,i) )  &
1062                                          * momentumflux_output_conversion(k)  &
1063                                                    * ust * rmask(j,i,sr)
1064!
1065!--               Momentum flux w*v*
1066                  sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5_wp *                 &
1067                                                    ( w(k,j-1,i) + w(k,j,i) )  &
1068                                          * momentumflux_output_conversion(k)  &
1069                                                    * vst * rmask(j,i,sr)
1070               ENDDO
1071            ENDDO
1072         ENDDO
1073
1074       ENDIF
1075       IF ( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
1076         !$OMP DO
1077         DO  i = nxl, nxr
1078            DO  j = nys, nyn
1079               DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
1080!
1081!--               Vertical heat flux
1082                  sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5_wp *                 &
1083                           ( pt(k,j,i)   - hom(k,1,4,sr) +                     &
1084                             pt(k+1,j,i) - hom(k+1,1,4,sr) )                   &
1085                           * heatflux_output_conversion(k)                     &
1086                           * w(k,j,i) * rmask(j,i,sr)
1087                  IF ( humidity )  THEN
1088                     pts = 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +            &
1089                                      q(k+1,j,i) - hom(k+1,1,41,sr) )
1090                     sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *      &
1091                                       waterflux_output_conversion(k) *        &
1092                                       rmask(j,i,sr)
1093                  ENDIF
1094                  IF ( passive_scalar )  THEN
1095                     pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,117,sr) +            &
1096                                      s(k+1,j,i) - hom(k+1,1,117,sr) )
1097                     sums_l(k,116,tn) = sums_l(k,116,tn) + pts * w(k,j,i) *     &
1098                                        rmask(j,i,sr)
1099                  ENDIF
1100               ENDDO
1101            ENDDO
1102         ENDDO
1103
1104       ENDIF
1105
1106!
1107!--    Density at top follows Neumann condition
1108       IF ( ocean )  THEN
1109          sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
1110          sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
1111       ENDIF
1112
1113!
1114!--    Divergence of vertical flux of resolved scale energy and pressure
1115!--    fluctuations as well as flux of pressure fluctuation itself (68).
1116!--    First calculate the products, then the divergence.
1117!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
1118       IF ( hom(nzb+1,2,55,0) /= 0.0_wp  .OR.  hom(nzb+1,2,68,0) /= 0.0_wp )   &
1119       THEN
1120          sums_ll = 0.0_wp  ! local array
1121
1122          !$OMP DO
1123          DO  i = nxl, nxr
1124             DO  j = nys, nyn
1125                DO  k = nzb_s_inner(j,i)+1, nzt
1126
1127                   sums_ll(k,1) = sums_ll(k,1) + 0.5_wp * w(k,j,i) * (         &
1128                  ( 0.25_wp * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1) )  &
1129                            - 0.5_wp * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) ) )**2&
1130                + ( 0.25_wp * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i) )  &
1131                            - 0.5_wp * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) ) )**2&
1132                + w(k,j,i)**2                                        )
1133
1134                   sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i)             &
1135                                               * ( p(k,j,i) + p(k+1,j,i) )
1136
1137                ENDDO
1138             ENDDO
1139          ENDDO
1140          sums_ll(0,1)     = 0.0_wp    ! because w is zero at the bottom
1141          sums_ll(nzt+1,1) = 0.0_wp
1142          sums_ll(0,2)     = 0.0_wp
1143          sums_ll(nzt+1,2) = 0.0_wp
1144
1145          DO  k = nzb+1, nzt
1146             sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
1147             sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
1148             sums_l(k,68,tn) = sums_ll(k,2)
1149          ENDDO
1150          sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
1151          sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
1152          sums_l(nzb,68,tn) = 0.0_wp    ! because w* = 0 at nzb
1153
1154       ENDIF
1155
1156!
1157!--    Divergence of vertical flux of SGS TKE and the flux itself (69)
1158       IF ( hom(nzb+1,2,57,0) /= 0.0_wp  .OR.  hom(nzb+1,2,69,0) /= 0.0_wp )   &
1159       THEN
1160          !$OMP DO
1161          DO  i = nxl, nxr
1162             DO  j = nys, nyn
1163                DO  k = nzb_s_inner(j,i)+1, nzt
1164
1165                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5_wp * (              &
1166                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
1167                 - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
1168                                                                ) * ddzw(k)
1169
1170                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5_wp * (              &
1171                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
1172                                                                )
1173
1174                ENDDO
1175             ENDDO
1176          ENDDO
1177          sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
1178          sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
1179
1180       ENDIF
1181
1182!
1183!--    Horizontal heat fluxes (subgrid, resolved, total).
1184!--    Do it only, if profiles shall be plotted.
1185       IF ( hom(nzb+1,2,58,0) /= 0.0_wp ) THEN
1186
1187          !$OMP DO
1188          DO  i = nxl, nxr
1189             DO  j = nys, nyn
1190                DO  k = nzb_s_inner(j,i)+1, nzt
1191!
1192!--                Subgrid horizontal heat fluxes u"pt", v"pt"
1193                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5_wp *                &
1194                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
1195                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
1196                                               * rho_air_zw(k)                 &
1197                                               * heatflux_output_conversion(k) &
1198                                                 * ddx * rmask(j,i,sr)
1199                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
1200                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
1201                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
1202                                               * rho_air_zw(k)                 &
1203                                               * heatflux_output_conversion(k) &
1204                                                 * ddy * rmask(j,i,sr)
1205!
1206!--                Resolved horizontal heat fluxes u*pt*, v*pt*
1207                   sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
1208                                                  ( u(k,j,i) - hom(k,1,1,sr) ) &
1209                                    * 0.5_wp * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
1210                                                 pt(k,j,i)   - hom(k,1,4,sr) ) &
1211                                               * heatflux_output_conversion(k)
1212                   pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
1213                                    pt(k,j,i)   - hom(k,1,4,sr) )
1214                   sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
1215                                                  ( v(k,j,i) - hom(k,1,2,sr) ) &
1216                                    * 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
1217                                                 pt(k,j,i)   - hom(k,1,4,sr) ) &
1218                                               * heatflux_output_conversion(k)
1219                ENDDO
1220             ENDDO
1221          ENDDO
1222!
1223!--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
1224          sums_l(nzb,58,tn) = 0.0_wp
1225          sums_l(nzb,59,tn) = 0.0_wp
1226          sums_l(nzb,60,tn) = 0.0_wp
1227          sums_l(nzb,61,tn) = 0.0_wp
1228          sums_l(nzb,62,tn) = 0.0_wp
1229          sums_l(nzb,63,tn) = 0.0_wp
1230
1231       ENDIF
1232       !$OMP END PARALLEL
1233
1234!
1235!--    Collect current large scale advection and subsidence tendencies for
1236!--    data output
1237       IF ( large_scale_forcing  .AND.  ( simulated_time > 0.0_wp ) )  THEN
1238!
1239!--       Interpolation in time of LSF_DATA
1240          nt = 1
1241          DO WHILE ( simulated_time - dt_3d > time_vert(nt) )
1242             nt = nt + 1
1243          ENDDO
1244          IF ( simulated_time - dt_3d /= time_vert(nt) )  THEN
1245            nt = nt - 1
1246          ENDIF
1247
1248          fac = ( simulated_time - dt_3d - time_vert(nt) )                     &
1249                / ( time_vert(nt+1)-time_vert(nt) )
1250
1251
1252          DO  k = nzb, nzt
1253             sums_ls_l(k,0) = td_lsa_lpt(k,nt)                                 &
1254                              + fac * ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )
1255             sums_ls_l(k,1) = td_lsa_q(k,nt)                                   &
1256                              + fac * ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) )
1257          ENDDO
1258
1259          sums_ls_l(nzt+1,0) = sums_ls_l(nzt,0)
1260          sums_ls_l(nzt+1,1) = sums_ls_l(nzt,1)
1261
1262          IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
1263
1264             DO  k = nzb, nzt
1265                sums_ls_l(k,2) = td_sub_lpt(k,nt) + fac *                      &
1266                                 ( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )
1267                sums_ls_l(k,3) = td_sub_q(k,nt) + fac *                        &
1268                                 ( td_sub_q(k,nt+1) - td_sub_q(k,nt) )
1269             ENDDO
1270
1271             sums_ls_l(nzt+1,2) = sums_ls_l(nzt,2)
1272             sums_ls_l(nzt+1,3) = sums_ls_l(nzt,3)
1273
1274          ENDIF
1275
1276       ENDIF
1277
1278       !$OMP PARALLEL PRIVATE( i, j, k, tn )
1279!$     tn = omp_get_thread_num()
1280       IF ( land_surface )  THEN
1281          !$OMP DO
1282          DO  i = nxl, nxr
1283             DO  j =  nys, nyn
1284                DO  k = nzb_soil, nzt_soil
1285                   sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil(k,j,i)         &
1286                                      * rmask(j,i,sr)
1287                   sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil(k,j,i)         &
1288                                      * rmask(j,i,sr)
1289                ENDDO
1290             ENDDO
1291          ENDDO
1292       ENDIF
1293       
1294       IF ( radiation .AND. radiation_scheme == 'rrtmg' )  THEN
1295          !$OMP DO
1296          DO  i = nxl, nxr
1297             DO  j =  nys, nyn
1298                DO  k = nzb_s_inner(j,i)+1, nzt+1
1299                   sums_l(k,102,tn)  = sums_l(k,102,tn)  + rad_lw_in(k,j,i)    &
1300                                       * rmask(j,i,sr)
1301                   sums_l(k,103,tn)  = sums_l(k,103,tn)  + rad_lw_out(k,j,i)   &
1302                                       * rmask(j,i,sr)
1303                   sums_l(k,104,tn)  = sums_l(k,104,tn)  + rad_sw_in(k,j,i)    &
1304                                       * rmask(j,i,sr)
1305                   sums_l(k,105,tn)  = sums_l(k,105,tn)  + rad_sw_out(k,j,i)   &
1306                                       * rmask(j,i,sr)
1307                   sums_l(k,106,tn)  = sums_l(k,106,tn)  + rad_lw_cs_hr(k,j,i) &
1308                                       * rmask(j,i,sr)
1309                   sums_l(k,107,tn)  = sums_l(k,107,tn)  + rad_lw_hr(k,j,i)    &
1310                                       * rmask(j,i,sr)
1311                   sums_l(k,108,tn)  = sums_l(k,108,tn)  + rad_sw_cs_hr(k,j,i) &
1312                                       * rmask(j,i,sr)
1313                   sums_l(k,109,tn)  = sums_l(k,109,tn)  + rad_sw_hr(k,j,i)    &
1314                                       * rmask(j,i,sr)
1315                ENDDO
1316             ENDDO
1317          ENDDO
1318       ENDIF
1319!
1320!--    Calculate the user-defined profiles
1321       CALL user_statistics( 'profiles', sr, tn )
1322       !$OMP END PARALLEL
1323
1324!
1325!--    Summation of thread sums
1326       IF ( threads_per_task > 1 )  THEN
1327          DO  i = 1, threads_per_task-1
1328             sums_l(:,3,0)          = sums_l(:,3,0) + sums_l(:,3,i)
1329             sums_l(:,4:40,0)       = sums_l(:,4:40,0) + sums_l(:,4:40,i)
1330             sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
1331                                      sums_l(:,45:pr_palm,i)
1332             IF ( max_pr_user > 0 )  THEN
1333                sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
1334                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
1335                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
1336             ENDIF
1337          ENDDO
1338       ENDIF
1339
1340#if defined( __parallel )
1341
1342!
1343!--    Compute total sum from local sums
1344       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1345       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL,   &
1346                           MPI_SUM, comm2d, ierr )
1347       IF ( large_scale_forcing )  THEN
1348          CALL MPI_ALLREDUCE( sums_ls_l(nzb,2), sums(nzb,83), ngp_sums_ls,     &
1349                              MPI_REAL, MPI_SUM, comm2d, ierr )
1350       ENDIF
1351#else
1352       sums = sums_l(:,:,0)
1353       IF ( large_scale_forcing )  THEN
1354          sums(:,81:88) = sums_ls_l
1355       ENDIF
1356#endif
1357
1358!
1359!--    Final values are obtained by division by the total number of grid points
1360!--    used for summation. After that store profiles.
1361!--    Check, if statistical regions do contain at least one grid point at the
1362!--    respective k-level, otherwise division by zero will lead to undefined
1363!--    values, which may cause e.g. problems with NetCDF output
1364!--    Profiles:
1365       DO  k = nzb, nzt+1
1366          sums(k,3)             = sums(k,3)             / ngp_2dh(sr)
1367          sums(k,12:22)         = sums(k,12:22)         / ngp_2dh(sr)
1368          sums(k,30:32)         = sums(k,30:32)         / ngp_2dh(sr)
1369          sums(k,35:39)         = sums(k,35:39)         / ngp_2dh(sr)
1370          sums(k,45:53)         = sums(k,45:53)         / ngp_2dh(sr)
1371          sums(k,55:63)         = sums(k,55:63)         / ngp_2dh(sr)
1372          sums(k,81:88)         = sums(k,81:88)         / ngp_2dh(sr)
1373          sums(k,89:114)        = sums(k,89:114)        / ngp_2dh(sr)
1374          sums(k,116)           = sums(k,116)           / ngp_2dh(sr)
1375          sums(k,119)           = sums(k,119)           / ngp_2dh(sr)
1376          IF ( ngp_2dh_s_inner(k,sr) /= 0 )  THEN
1377             sums(k,8:11)          = sums(k,8:11)          / ngp_2dh_s_inner(k,sr)
1378             sums(k,23:29)         = sums(k,23:29)         / ngp_2dh_s_inner(k,sr)
1379             sums(k,33:34)         = sums(k,33:34)         / ngp_2dh_s_inner(k,sr)
1380             sums(k,40)            = sums(k,40)            / ngp_2dh_s_inner(k,sr)
1381             sums(k,54)            = sums(k,54)            / ngp_2dh_s_inner(k,sr)
1382             sums(k,64)            = sums(k,64)            / ngp_2dh_s_inner(k,sr)
1383             sums(k,70:80)         = sums(k,70:80)         / ngp_2dh_s_inner(k,sr)
1384             sums(k,118)           = sums(k,118)           / ngp_2dh_s_inner(k,sr)
1385             sums(k,120:pr_palm-2) = sums(k,120:pr_palm-2) / ngp_2dh_s_inner(k,sr)
1386          ENDIF
1387       ENDDO
1388
1389!--    u* and so on
1390!--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
1391!--    size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
1392!--    above the topography, they are being divided by ngp_2dh(sr)
1393       sums(nzb:nzb+3,pr_palm)    = sums(nzb:nzb+3,pr_palm)    / &
1394                                    ngp_2dh(sr)
1395       sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)       / &    ! qs
1396                                    ngp_2dh(sr)
1397       sums(nzb+13,pr_palm)       = sums(nzb+13,pr_palm)       / &    ! ss
1398                                    ngp_2dh(sr)
1399!--    eges, e*
1400       sums(nzb+4:nzb+5,pr_palm)  = sums(nzb+4:nzb+5,pr_palm)  / &
1401                                    ngp_3d(sr)
1402!--    Old and new divergence
1403       sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
1404                                    ngp_3d_inner(sr)
1405
1406!--    User-defined profiles
1407       IF ( max_pr_user > 0 )  THEN
1408          DO  k = nzb, nzt+1
1409             sums(k,pr_palm+1:pr_palm+max_pr_user) = &
1410                                    sums(k,pr_palm+1:pr_palm+max_pr_user) / &
1411                                    ngp_2dh_s_inner(k,sr)
1412          ENDDO
1413       ENDIF
1414
1415!
1416!--    Collect horizontal average in hom.
1417!--    Compute deduced averages (e.g. total heat flux)
1418       hom(:,1,3,sr)  = sums(:,3)      ! w
1419       hom(:,1,8,sr)  = sums(:,8)      ! e     profiles 5-7 are initial profiles
1420       hom(:,1,9,sr)  = sums(:,9)      ! km
1421       hom(:,1,10,sr) = sums(:,10)     ! kh
1422       hom(:,1,11,sr) = sums(:,11)     ! l
1423       hom(:,1,12,sr) = sums(:,12)     ! w"u"
1424       hom(:,1,13,sr) = sums(:,13)     ! w*u*
1425       hom(:,1,14,sr) = sums(:,14)     ! w"v"
1426       hom(:,1,15,sr) = sums(:,15)     ! w*v*
1427       hom(:,1,16,sr) = sums(:,16)     ! w"pt"
1428       hom(:,1,17,sr) = sums(:,17)     ! w*pt*
1429       hom(:,1,18,sr) = sums(:,16) + sums(:,17)    ! wpt
1430       hom(:,1,19,sr) = sums(:,12) + sums(:,13)    ! wu
1431       hom(:,1,20,sr) = sums(:,14) + sums(:,15)    ! wv
1432       hom(:,1,21,sr) = sums(:,21)     ! w*pt*BC
1433       hom(:,1,22,sr) = sums(:,16) + sums(:,21)    ! wptBC
1434                                       ! profile 24 is initial profile (sa)
1435                                       ! profiles 25-29 left empty for initial
1436                                       ! profiles
1437       hom(:,1,30,sr) = sums(:,30)     ! u*2
1438       hom(:,1,31,sr) = sums(:,31)     ! v*2
1439       hom(:,1,32,sr) = sums(:,32)     ! w*2
1440       hom(:,1,33,sr) = sums(:,33)     ! pt*2
1441       hom(:,1,34,sr) = sums(:,34)     ! e*
1442       hom(:,1,35,sr) = sums(:,35)     ! w*2pt*
1443       hom(:,1,36,sr) = sums(:,36)     ! w*pt*2
1444       hom(:,1,37,sr) = sums(:,37)     ! w*e*
1445       hom(:,1,38,sr) = sums(:,38)     ! w*3
1446       hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20_wp )**1.5_wp   ! Sw
1447       hom(:,1,40,sr) = sums(:,40)     ! p
1448       hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
1449       hom(:,1,46,sr) = sums(:,46)     ! w*vpt*       
1450       hom(:,1,47,sr) = sums(:,45) + sums(:,46)    ! wvpt
1451       hom(:,1,48,sr) = sums(:,48)     ! w"q" (w"qv")
1452       hom(:,1,49,sr) = sums(:,49)     ! w*q* (w*qv*)
1453       hom(:,1,50,sr) = sums(:,48) + sums(:,49)    ! wq (wqv)
1454       hom(:,1,51,sr) = sums(:,51)     ! w"qv"
1455       hom(:,1,52,sr) = sums(:,52)     ! w*qv*       
1456       hom(:,1,53,sr) = sums(:,52) + sums(:,51)    ! wq (wqv)
1457       hom(:,1,54,sr) = sums(:,54)     ! ql
1458       hom(:,1,55,sr) = sums(:,55)     ! w*u*u*/dz
1459       hom(:,1,56,sr) = sums(:,56)     ! w*p*/dz
1460       hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho_ocean )/dz
1461       hom(:,1,58,sr) = sums(:,58)     ! u"pt"
1462       hom(:,1,59,sr) = sums(:,59)     ! u*pt*
1463       hom(:,1,60,sr) = sums(:,58) + sums(:,59)    ! upt_t
1464       hom(:,1,61,sr) = sums(:,61)     ! v"pt"
1465       hom(:,1,62,sr) = sums(:,62)     ! v*pt*
1466       hom(:,1,63,sr) = sums(:,61) + sums(:,62)    ! vpt_t
1467       hom(:,1,64,sr) = sums(:,64)     ! rho_ocean
1468       hom(:,1,65,sr) = sums(:,65)     ! w"sa"
1469       hom(:,1,66,sr) = sums(:,66)     ! w*sa*
1470       hom(:,1,67,sr) = sums(:,65) + sums(:,66)    ! wsa
1471       hom(:,1,68,sr) = sums(:,68)     ! w*p*
1472       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho_ocean
1473       hom(:,1,70,sr) = sums(:,70)     ! q*2
1474       hom(:,1,71,sr) = sums(:,71)     ! prho
1475       hom(:,1,72,sr) = hyp * 1E-4_wp  ! hyp in dbar
1476       hom(:,1,73,sr) = sums(:,73)     ! nr
1477       hom(:,1,74,sr) = sums(:,74)     ! qr
1478       hom(:,1,75,sr) = sums(:,75)     ! qc
1479       hom(:,1,76,sr) = sums(:,76)     ! prr (precipitation rate)
1480                                       ! 77 is initial density profile
1481       hom(:,1,78,sr) = ug             ! ug
1482       hom(:,1,79,sr) = vg             ! vg
1483       hom(:,1,80,sr) = w_subs         ! w_subs
1484
1485       IF ( large_scale_forcing )  THEN
1486          hom(:,1,81,sr) = sums_ls_l(:,0)          ! td_lsa_lpt
1487          hom(:,1,82,sr) = sums_ls_l(:,1)          ! td_lsa_q
1488          IF ( use_subsidence_tendencies )  THEN
1489             hom(:,1,83,sr) = sums_ls_l(:,2)       ! td_sub_lpt
1490             hom(:,1,84,sr) = sums_ls_l(:,3)       ! td_sub_q
1491          ELSE
1492             hom(:,1,83,sr) = sums(:,83)           ! td_sub_lpt
1493             hom(:,1,84,sr) = sums(:,84)           ! td_sub_q
1494          ENDIF
1495          hom(:,1,85,sr) = sums(:,85)              ! td_nud_lpt
1496          hom(:,1,86,sr) = sums(:,86)              ! td_nud_q
1497          hom(:,1,87,sr) = sums(:,87)              ! td_nud_u
1498          hom(:,1,88,sr) = sums(:,88)              ! td_nud_v
1499       ENDIF
1500
1501       IF ( land_surface )  THEN
1502          hom(:,1,89,sr) = sums(:,89)              ! t_soil
1503                                                   ! 90 is initial t_soil profile
1504          hom(:,1,91,sr) = sums(:,91)              ! m_soil
1505                                                   ! 92 is initial m_soil profile
1506          hom(:,1,93,sr)  = sums(:,93)             ! ghf_eb
1507          hom(:,1,94,sr)  = sums(:,94)             ! shf_eb
1508          hom(:,1,95,sr)  = sums(:,95)             ! qsws_eb
1509          hom(:,1,96,sr)  = sums(:,96)             ! qsws_liq_eb
1510          hom(:,1,97,sr)  = sums(:,97)             ! qsws_soil_eb
1511          hom(:,1,98,sr)  = sums(:,98)             ! qsws_veg_eb
1512          hom(:,1,99,sr)  = sums(:,99)             ! r_a
1513          hom(:,1,100,sr) = sums(:,100)            ! r_s
1514
1515       ENDIF
1516
1517       IF ( radiation )  THEN
1518          hom(:,1,101,sr) = sums(:,101)            ! rad_net
1519          hom(:,1,102,sr) = sums(:,102)            ! rad_lw_in
1520          hom(:,1,103,sr) = sums(:,103)            ! rad_lw_out
1521          hom(:,1,104,sr) = sums(:,104)            ! rad_sw_in
1522          hom(:,1,105,sr) = sums(:,105)            ! rad_sw_out
1523
1524          IF ( radiation_scheme == 'rrtmg' )  THEN
1525             hom(:,1,106,sr) = sums(:,106)            ! rad_lw_cs_hr
1526             hom(:,1,107,sr) = sums(:,107)            ! rad_lw_hr
1527             hom(:,1,108,sr) = sums(:,108)            ! rad_sw_cs_hr
1528             hom(:,1,109,sr) = sums(:,109)            ! rad_sw_hr
1529
1530             hom(:,1,110,sr) = sums(:,110)            ! rrtm_aldif
1531             hom(:,1,111,sr) = sums(:,111)            ! rrtm_aldir
1532             hom(:,1,112,sr) = sums(:,112)            ! rrtm_asdif
1533             hom(:,1,113,sr) = sums(:,113)            ! rrtm_asdir
1534          ENDIF
1535       ENDIF
1536
1537       hom(:,1,114,sr) = sums(:,114)            !: L
1538
1539       IF ( passive_scalar )  THEN
1540          hom(:,1,119,sr) = sums(:,119)     ! w"s"
1541          hom(:,1,116,sr) = sums(:,116)     ! w*s*
1542          hom(:,1,120,sr) = sums(:,119) + sums(:,116)    ! ws
1543          hom(:,1,118,sr) = sums(:,118)     ! s*2
1544       ENDIF
1545
1546       hom(:,1,121,sr) = rho_air       ! rho_air in Kg/m^3
1547       hom(:,1,122,sr) = rho_air_zw    ! rho_air_zw in Kg/m^3
1548
1549       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
1550                                       ! u*, w'u', w'v', t* (in last profile)
1551
1552       IF ( max_pr_user > 0 )  THEN    ! user-defined profiles
1553          hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
1554                               sums(:,pr_palm+1:pr_palm+max_pr_user)
1555       ENDIF
1556
1557!
1558!--    Determine the boundary layer height using two different schemes.
1559!--    First scheme: Starting from the Earth's (Ocean's) surface, look for the
1560!--    first relative minimum (maximum) of the total heat flux.
1561!--    The corresponding height is assumed as the boundary layer height, if it
1562!--    is less than 1.5 times the height where the heat flux becomes negative
1563!--    (positive) for the first time.
1564       z_i(1) = 0.0_wp
1565       first = .TRUE.
1566
1567       IF ( ocean )  THEN
1568          DO  k = nzt, nzb+1, -1
1569             IF ( first  .AND.  hom(k,1,18,sr) < -1.0E-8_wp )  THEN
1570                first = .FALSE.
1571                height = zw(k)
1572             ENDIF
1573             IF ( hom(k,1,18,sr) < -1.0E-8_wp  .AND.                           &
1574                  hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
1575                IF ( zw(k) < 1.5_wp * height )  THEN
1576                   z_i(1) = zw(k)
1577                ELSE
1578                   z_i(1) = height
1579                ENDIF
1580                EXIT
1581             ENDIF
1582          ENDDO
1583       ELSE
1584          DO  k = nzb, nzt-1
1585             IF ( first  .AND.  hom(k,1,18,sr) < -1.0E-8_wp )  THEN
1586                first = .FALSE.
1587                height = zw(k)
1588             ENDIF
1589             IF ( hom(k,1,18,sr) < -1.0E-8_wp  .AND.                           &
1590                  hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
1591                IF ( zw(k) < 1.5_wp * height )  THEN
1592                   z_i(1) = zw(k)
1593                ELSE
1594                   z_i(1) = height
1595                ENDIF
1596                EXIT
1597             ENDIF
1598          ENDDO
1599       ENDIF
1600
1601!
1602!--    Second scheme: Gradient scheme from Sullivan et al. (1998), modified
1603!--    by Uhlenbrock(2006). The boundary layer height is the height with the
1604!--    maximal local temperature gradient: starting from the second (the last
1605!--    but one) vertical gridpoint, the local gradient must be at least
1606!--    0.2K/100m and greater than the next four gradients.
1607!--    WARNING: The threshold value of 0.2K/100m must be adjusted for the
1608!--             ocean case!
1609       z_i(2) = 0.0_wp
1610       DO  k = nzb+1, nzt+1
1611          dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
1612       ENDDO
1613       dptdz_threshold = 0.2_wp / 100.0_wp
1614
1615       IF ( ocean )  THEN
1616          DO  k = nzt+1, nzb+5, -1
1617             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
1618                  dptdz(k) > dptdz(k-1)  .AND.  dptdz(k) > dptdz(k-2)  .AND.  &
1619                  dptdz(k) > dptdz(k-3)  .AND.  dptdz(k) > dptdz(k-4) )  THEN
1620                z_i(2) = zw(k-1)
1621                EXIT
1622             ENDIF
1623          ENDDO
1624       ELSE
1625          DO  k = nzb+1, nzt-3
1626             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
1627                  dptdz(k) > dptdz(k+1)  .AND.  dptdz(k) > dptdz(k+2)  .AND.  &
1628                  dptdz(k) > dptdz(k+3)  .AND.  dptdz(k) > dptdz(k+4) )  THEN
1629                z_i(2) = zw(k-1)
1630                EXIT
1631             ENDIF
1632          ENDDO
1633       ENDIF
1634
1635       hom(nzb+6,1,pr_palm,sr) = z_i(1)
1636       hom(nzb+7,1,pr_palm,sr) = z_i(2)
1637
1638!
1639!--    Determine vertical index which is nearest to the mean surface level
1640!--    height of the respective statistic region
1641       DO  k = nzb, nzt
1642          IF ( zw(k) >= mean_surface_level_height(sr) )  THEN
1643             k_surface_level = k
1644             EXIT
1645          ENDIF
1646       ENDDO
1647!
1648!--    Computation of both the characteristic vertical velocity and
1649!--    the characteristic convective boundary layer temperature.
1650!--    The inversion height entering into the equation is defined with respect
1651!--    to the mean surface level height of the respective statistic region.
1652!--    The horizontal average at surface level index + 1 is input for the
1653!--    average temperature.
1654       IF ( hom(k_surface_level,1,18,sr) > 1.0E-8_wp  .AND.  z_i(1) /= 0.0_wp )&
1655       THEN
1656          hom(nzb+8,1,pr_palm,sr) = &
1657             ( g / hom(k_surface_level+1,1,4,sr) *                             &
1658             ( hom(k_surface_level,1,18,sr) / heatflux_output_conversion(nzb) )&
1659             * ABS( z_i(1) - mean_surface_level_height(sr) ) )**0.333333333_wp
1660       ELSE
1661          hom(nzb+8,1,pr_palm,sr)  = 0.0_wp
1662       ENDIF
1663
1664!
1665!--    Collect the time series quantities
1666       ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)     ! E
1667       ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)     ! E*
1668       ts_value(3,sr) = dt_3d
1669       ts_value(4,sr) = hom(nzb,1,pr_palm,sr)       ! u*
1670       ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr)     ! th*
1671       ts_value(6,sr) = u_max
1672       ts_value(7,sr) = v_max
1673       ts_value(8,sr) = w_max
1674       ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr)    ! new divergence
1675       ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr)    ! old Divergence
1676       ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr)    ! z_i(1)
1677       ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr)    ! z_i(2)
1678       ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr)    ! w*
1679       ts_value(14,sr) = hom(nzb,1,16,sr)           ! w'pt'   at k=0
1680       ts_value(15,sr) = hom(nzb+1,1,16,sr)         ! w'pt'   at k=1
1681       ts_value(16,sr) = hom(nzb+1,1,18,sr)         ! wpt     at k=1
1682       ts_value(17,sr) = hom(nzb,1,4,sr)            ! pt(0)
1683       ts_value(18,sr) = hom(nzb+1,1,4,sr)          ! pt(zp)
1684       ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr)    ! u'w'    at k=0
1685       ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)    ! v'w'    at k=0
1686       ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
1687
1688       IF ( .NOT. neutral )  THEN
1689          ts_value(22,sr) = hom(nzb,1,114,sr)          ! L
1690       ELSE
1691          ts_value(22,sr) = 1.0E10_wp
1692       ENDIF
1693
1694       ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
1695
1696       IF ( passive_scalar )  THEN
1697          ts_value(24,sr) = hom(nzb+13,1,119,sr)       ! w"s" ( to do ! )
1698          ts_value(25,sr) = hom(nzb+13,1,pr_palm,sr)   ! s*
1699       ENDIF
1700
1701!
1702!--    Collect land surface model timeseries
1703       IF ( land_surface )  THEN
1704          ts_value(dots_soil  ,sr) = hom(nzb,1,93,sr)           ! ghf_eb
1705          ts_value(dots_soil+1,sr) = hom(nzb,1,94,sr)           ! shf_eb
1706          ts_value(dots_soil+2,sr) = hom(nzb,1,95,sr)           ! qsws_eb
1707          ts_value(dots_soil+3,sr) = hom(nzb,1,96,sr)           ! qsws_liq_eb
1708          ts_value(dots_soil+4,sr) = hom(nzb,1,97,sr)           ! qsws_soil_eb
1709          ts_value(dots_soil+5,sr) = hom(nzb,1,98,sr)           ! qsws_veg_eb
1710          ts_value(dots_soil+6,sr) = hom(nzb,1,99,sr)           ! r_a
1711          ts_value(dots_soil+7,sr) = hom(nzb,1,100,sr)          ! r_s
1712       ENDIF
1713!
1714!--    Collect radiation model timeseries
1715       IF ( radiation )  THEN
1716          ts_value(dots_rad,sr)   = hom(nzb,1,101,sr)          ! rad_net
1717          ts_value(dots_rad+1,sr) = hom(nzb,1,102,sr)          ! rad_lw_in
1718          ts_value(dots_rad+2,sr) = hom(nzb,1,103,sr)          ! rad_lw_out
1719          ts_value(dots_rad+3,sr) = hom(nzb,1,104,sr)          ! rad_sw_in
1720          ts_value(dots_rad+4,sr) = hom(nzb,1,105,sr)          ! rad_sw_out
1721
1722          IF ( radiation_scheme == 'rrtmg' )  THEN
1723             ts_value(dots_rad+5,sr) = hom(nzb,1,110,sr)          ! rrtm_aldif
1724             ts_value(dots_rad+6,sr) = hom(nzb,1,111,sr)          ! rrtm_aldir
1725             ts_value(dots_rad+7,sr) = hom(nzb,1,112,sr)          ! rrtm_asdif
1726             ts_value(dots_rad+8,sr) = hom(nzb,1,113,sr)          ! rrtm_asdir
1727          ENDIF
1728
1729       ENDIF
1730
1731!
1732!--    Calculate additional statistics provided by the user interface
1733       CALL user_statistics( 'time_series', sr, 0 )
1734
1735    ENDDO    ! loop of the subregions
1736
1737!
1738!-- If required, sum up horizontal averages for subsequent time averaging.
1739!-- Do not sum, if flow statistics is called before the first initial time step.
1740    IF ( do_sum  .AND.  simulated_time /= 0.0_wp )  THEN
1741       IF ( average_count_pr == 0 )  hom_sum = 0.0_wp
1742       hom_sum = hom_sum + hom(:,1,:,:)
1743       average_count_pr = average_count_pr + 1
1744       do_sum = .FALSE.
1745    ENDIF
1746
1747!
1748!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
1749!-- This flag is reset after each time step in time_integration.
1750    flow_statistics_called = .TRUE.
1751
1752    CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
1753
1754
1755 END SUBROUTINE flow_statistics
Note: See TracBrowser for help on using the repository browser.