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

Last change on this file since 2512 was 2320, checked in by suehring, 7 years ago

large-scale forcing and nudging modularized

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