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

Last change on this file since 2258 was 2252, checked in by knoop, 7 years ago

air density now depending on surface_pressure even in boussinesq mode

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