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

Last change on this file since 2186 was 2119, checked in by raasch, 8 years ago

last commit documented

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