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

Last change on this file since 3003 was 3003, checked in by Giersch, 3 years ago

Bugfix: W* and Z_I in the first line of the run control file of restarts correspond now to the values in the last run control output line of the previous run

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