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

Last change on this file since 2716 was 2716, checked in by kanani, 6 years ago

Correction of "Former revisions" section

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