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

Last change on this file since 2292 was 2292, checked in by schwenkel, 4 years ago

implementation of new bulk microphysics scheme

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