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

Last change on this file since 2296 was 2296, checked in by maronga, 4 years ago

added new spinup mechanism for surface/radiation models

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