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

Last change on this file since 2773 was 2773, checked in by suehring, 3 years ago

Nesting for chemical species implemented; Bugfix passive scalar boundary condition after anterpolation; Timeseries output of surface temperature; Enable initialization of 3D topography (was commented out so far)

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