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

Last change on this file since 2699 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk

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