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

Last change on this file since 3046 was 3042, checked in by schwenkel, 3 years ago

Changed the name specific humidity to mixing ratio

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