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

Last change on this file since 1106 was 1054, checked in by hoffmann, 12 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 56.2 KB
Line 
1 SUBROUTINE flow_statistics
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later 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-2012  Leibniz University Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23! Former revisions:
24! -----------------
25! $Id: flow_statistics.f90 1054 2012-11-13 17:30:09Z raasch $
26!
27! 1053 2012-11-13 17:11:03Z hoffmann
28! necessary additions for two-moment cloud physics scheme:
29! +nr, qr, qc, prr
30!
31! 1036 2012-10-22 13:43:42Z raasch
32! code put under GPL (PALM 3.9)
33!
34! 1007 2012-09-19 14:30:36Z franke
35! Calculation of buoyancy flux for humidity in case of WS-scheme is now using
36! turbulent fluxes of WS-scheme
37! Bugfix: Calculation of subgridscale buoyancy flux for humidity and cloud
38! droplets at nzb and nzt added
39!
40! 801 2012-01-10 17:30:36Z suehring
41! Calculation of turbulent fluxes in advec_ws is now thread-safe.
42!
43! 743 2011-08-18 16:10:16Z suehring
44! Calculation of turbulent fluxes with WS-scheme only for the whole model
45! domain, not for user-defined subregions.
46!
47! 709 2011-03-30 09:31:40Z raasch
48! formatting adjustments
49!
50! 699 2011-03-22 17:52:22Z suehring
51! Bugfix in calculation of vertical velocity skewness. The added absolute value
52! avoid negative values in the root. Negative values of w'w' can occur at the
53! top or bottom of the model domain due to degrading the order of advection
54! scheme. Furthermore the calculation will be the same for all advection
55! schemes.
56!
57! 696 2011-03-18 07:03:49Z raasch
58! Bugfix: Summation of Wicker-Skamarock scheme fluxes and variances for all
59! threads
60!
61! 678 2011-02-02 14:31:56Z raasch
62! Bugfix in calculation of the divergence of vertical flux of resolved scale
63! energy, pressure fluctuations, and flux of pressure fluctuation itself
64!
65! 673 2011-01-18 16:19:48Z suehring
66! Top bc for the horizontal velocity variances added for ocean runs.
67! Setting the corresponding bottom bc moved to advec_ws.
68!
69! 667 2010-12-23 12:06:00Z suehring/gryschka
70! When advection is computed with ws-scheme, turbulent fluxes are already
71! computed in the respective advection routines and buffered in arrays
72! sums_xx_ws_l(). This is due to a consistent treatment of statistics with the
73! numerics and to avoid unphysical kinks near the surface.
74! So some if requests has to be done to dicern between fluxes from ws-scheme
75! other advection schemes.
76! Furthermore the computation of z_i is only done if the heat flux exceeds a
77! minimum value. This affects only simulations of a neutral boundary layer and
78! is due to reasons of computations in the advection scheme.
79!
80! 624 2010-12-10 11:46:30Z heinze
81! Calculation of q*2 added
82!
83! 622 2010-12-10 08:08:13Z raasch
84! optional barriers included in order to speed up collective operations
85!
86! 388 2009-09-23 09:40:33Z raasch
87! Vertical profiles of potential density and hydrostatic pressure are
88! calculated.
89! Added missing timeseries calculation of w"q"(0), moved timeseries q* to the
90! end.
91! Temperature gradient criterion for estimating the boundary layer height
92! replaced by the gradient criterion of Sullivan et al. (1998).
93! Output of messages replaced by message handling routine.
94!
95! 197 2008-09-16 15:29:03Z raasch
96! Spline timeseries splptx etc. removed, timeseries w'u', w'v', w'q' (k=0)
97! added,
98! bugfix: divide sums(k,8) (e) and sums(k,34) (e*) by ngp_2dh_s_inner(k,sr)
99! (like other scalars)
100!
101! 133 2007-11-20 10:10:53Z letzel
102! Vertical profiles now based on nzb_s_inner; they are divided by
103! ngp_2dh_s_inner (scalars, procucts of scalars) and ngp_2dh (staggered
104! velocity components and their products, procucts of scalars and velocity
105! components), respectively.
106!
107! 106 2007-08-16 14:30:26Z raasch
108! Prescribed momentum fluxes at the top surface are used,
109! profiles for w*p* and w"e are calculated
110!
111! 97 2007-06-21 08:23:15Z raasch
112! Statistics for ocean version (salinity, density) added,
113! calculation of z_i and Deardorff velocity scale adjusted to be used with
114! the ocean version
115!
116! 87 2007-05-22 15:46:47Z raasch
117! Two more arguments added to user_statistics, which is now also called for
118! user-defined profiles,
119! var_hom and var_sum renamed pr_palm
120!
121! 82 2007-04-16 15:40:52Z raasch
122! Cpp-directive lcmuk changed to intel_openmp_bug
123!
124! 75 2007-03-22 09:54:05Z raasch
125! Collection of time series quantities moved from routine flow_statistics to
126! here, routine user_statistics is called for each statistic region,
127! moisture renamed humidity
128!
129! 19 2007-02-23 04:53:48Z raasch
130! fluxes at top modified (tswst, qswst)
131!
132! RCS Log replace by Id keyword, revision history cleaned up
133!
134! Revision 1.41  2006/08/04 14:37:50  raasch
135! Error removed in non-parallel part (sums_l)
136!
137! Revision 1.1  1997/08/11 06:15:17  raasch
138! Initial revision
139!
140!
141! Description:
142! ------------
143! Compute average profiles and further average flow quantities for the different
144! user-defined (sub-)regions. The region indexed 0 is the total model domain.
145!
146! NOTE: For simplicity, nzb_s_inner and nzb_diff_s_inner are being used as a
147! ----  lower vertical index for k-loops for all variables, although strictly
148! speaking the k-loops would have to be split up according to the staggered
149! grid. However, this implies no error since staggered velocity components are
150! zero at the walls and inside buildings.
151!------------------------------------------------------------------------------!
152
153    USE arrays_3d
154    USE cloud_parameters
155    USE control_parameters
156    USE cpulog
157    USE grid_variables
158    USE indices
159    USE interfaces
160    USE pegrid
161    USE statistics
162
163    IMPLICIT NONE
164
165    INTEGER ::  i, j, k, omp_get_thread_num, sr, tn
166    LOGICAL ::  first
167    REAL    ::  dptdz_threshold, height, pts, sums_l_eper, sums_l_etot, ust, &
168                ust2, u2, vst, vst2, v2, w2, z_i(2)
169    REAL    ::  dptdz(nzb+1:nzt+1)
170    REAL    ::  sums_ll(nzb:nzt+1,2)
171
172    CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
173
174!
175!-- To be on the safe side, check whether flow_statistics has already been
176!-- called once after the current time step
177    IF ( flow_statistics_called )  THEN
178
179       message_string = 'flow_statistics is called two times within one ' // &
180                        'timestep'
181       CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
182
183    ENDIF
184
185!
186!-- Compute statistics for each (sub-)region
187    DO  sr = 0, statistic_regions
188
189!
190!--    Initialize (local) summation array
191       sums_l = 0.0
192
193!
194!--    Store sums that have been computed in other subroutines in summation
195!--    array
196       sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
197!--    WARNING: next line still has to be adjusted for OpenMP
198       sums_l(:,21,0) = sums_wsts_bc_l(:,sr)  ! heat flux from advec_s_bc
199       sums_l(nzb+9,pr_palm,0)  = sums_divold_l(sr)  ! old divergence from pres
200       sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
201
202!
203!--    Copy the turbulent quantities, evaluated in the advection routines to
204!--    the local array sums_l() for further computations
205       IF ( ws_scheme_mom .AND. sr == 0 )  THEN
206
207!
208!--       According to the Neumann bc for the horizontal velocity components,
209!--       the corresponding fluxes has to satisfiy the same bc.
210          IF ( ocean )  THEN
211             sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
212             sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)
213          ENDIF
214
215          DO  i = 0, threads_per_task-1
216!
217!--          Swap the turbulent quantities evaluated in advec_ws.
218             sums_l(:,13,i) = sums_wsus_ws_l(:,i)       ! w*u*
219             sums_l(:,15,i) = sums_wsvs_ws_l(:,i)       ! w*v*
220             sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2
221             sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2
222             sums_l(:,32,i) = sums_ws2_ws_l(:,i)        ! w*2
223             sums_l(:,34,i) = sums_l(:,34,i) + 0.5 *                        & 
224                              ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +   &
225                                sums_ws2_ws_l(:,i) )    ! e*
226             DO  k = nzb, nzt
227                sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5 * ( &
228                                                      sums_us2_ws_l(k,i) +  &
229                                                      sums_vs2_ws_l(k,i) +  &
230                                                      sums_ws2_ws_l(k,i) )
231             ENDDO
232          ENDDO
233
234       ENDIF
235
236       IF ( ws_scheme_sca .AND. sr == 0 )  THEN
237
238          DO  i = 0, threads_per_task-1
239             sums_l(:,17,i) = sums_wspts_ws_l(:,i)      ! w*pt* from advec_s_ws
240             IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa*
241             IF ( humidity .OR. passive_scalar ) sums_l(:,49,i) =              &
242                                                   sums_wsqs_ws_l(:,i) !w*q*
243          ENDDO
244
245       ENDIF
246!
247!--    Horizontally averaged profiles of horizontal velocities and temperature.
248!--    They must have been computed before, because they are already required
249!--    for other horizontal averages.
250       tn = 0
251
252       !$OMP PARALLEL PRIVATE( i, j, k, tn )
253#if defined( __intel_openmp_bug )
254       tn = omp_get_thread_num()
255#else
256!$     tn = omp_get_thread_num()
257#endif
258
259       !$OMP DO
260       DO  i = nxl, nxr
261          DO  j =  nys, nyn
262             DO  k = nzb_s_inner(j,i), nzt+1
263                sums_l(k,1,tn)  = sums_l(k,1,tn)  + u(k,j,i)  * rmask(j,i,sr)
264                sums_l(k,2,tn)  = sums_l(k,2,tn)  + v(k,j,i)  * rmask(j,i,sr)
265                sums_l(k,4,tn)  = sums_l(k,4,tn)  + pt(k,j,i) * rmask(j,i,sr)
266             ENDDO
267          ENDDO
268       ENDDO
269
270!
271!--    Horizontally averaged profile of salinity
272       IF ( ocean )  THEN
273          !$OMP DO
274          DO  i = nxl, nxr
275             DO  j =  nys, nyn
276                DO  k = nzb_s_inner(j,i), nzt+1
277                   sums_l(k,23,tn)  = sums_l(k,23,tn) + &
278                                      sa(k,j,i) * rmask(j,i,sr)
279                ENDDO
280             ENDDO
281          ENDDO
282       ENDIF
283
284!
285!--    Horizontally averaged profiles of virtual potential temperature,
286!--    total water content, specific humidity and liquid water potential
287!--    temperature
288       IF ( humidity )  THEN
289          !$OMP DO
290          DO  i = nxl, nxr
291             DO  j =  nys, nyn
292                DO  k = nzb_s_inner(j,i), nzt+1
293                   sums_l(k,44,tn)  = sums_l(k,44,tn) + &
294                                      vpt(k,j,i) * rmask(j,i,sr)
295                   sums_l(k,41,tn)  = sums_l(k,41,tn) + &
296                                      q(k,j,i) * rmask(j,i,sr)
297                ENDDO
298             ENDDO
299          ENDDO
300          IF ( cloud_physics )  THEN
301             !$OMP DO
302             DO  i = nxl, nxr
303                DO  j =  nys, nyn
304                   DO  k = nzb_s_inner(j,i), nzt+1
305                      sums_l(k,42,tn) = sums_l(k,42,tn) + &
306                                      ( q(k,j,i) - ql(k,j,i) ) * rmask(j,i,sr)
307                      sums_l(k,43,tn) = sums_l(k,43,tn) + ( &
308                                      pt(k,j,i) + l_d_cp*pt_d_t(k) * ql(k,j,i) &
309                                                          ) * rmask(j,i,sr)
310                   ENDDO
311                ENDDO
312             ENDDO
313          ENDIF
314       ENDIF
315
316!
317!--    Horizontally averaged profiles of passive scalar
318       IF ( passive_scalar )  THEN
319          !$OMP DO
320          DO  i = nxl, nxr
321             DO  j =  nys, nyn
322                DO  k = nzb_s_inner(j,i), nzt+1
323                   sums_l(k,41,tn)  = sums_l(k,41,tn) + q(k,j,i) * rmask(j,i,sr)
324                ENDDO
325             ENDDO
326          ENDDO
327       ENDIF
328       !$OMP END PARALLEL
329!
330!--    Summation of thread sums
331       IF ( threads_per_task > 1 )  THEN
332          DO  i = 1, threads_per_task-1
333             sums_l(:,1,0) = sums_l(:,1,0) + sums_l(:,1,i)
334             sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i)
335             sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i)
336             IF ( ocean )  THEN
337                sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)
338             ENDIF
339             IF ( humidity )  THEN
340                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
341                sums_l(:,44,0) = sums_l(:,44,0) + sums_l(:,44,i)
342                IF ( cloud_physics )  THEN
343                   sums_l(:,42,0) = sums_l(:,42,0) + sums_l(:,42,i)
344                   sums_l(:,43,0) = sums_l(:,43,0) + sums_l(:,43,i)
345                ENDIF
346             ENDIF
347             IF ( passive_scalar )  THEN
348                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
349             ENDIF
350          ENDDO
351       ENDIF
352
353#if defined( __parallel )
354!
355!--    Compute total sum from local sums
356       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
357       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL, &
358                           MPI_SUM, comm2d, ierr )
359       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
360       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL, &
361                           MPI_SUM, comm2d, ierr )
362       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
363       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL, &
364                           MPI_SUM, comm2d, ierr )
365       IF ( ocean )  THEN
366          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
367          CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb, &
368                              MPI_REAL, MPI_SUM, comm2d, ierr )
369       ENDIF
370       IF ( humidity ) THEN
371          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
372          CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb, &
373                              MPI_REAL, MPI_SUM, comm2d, ierr )
374          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
375          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
376                              MPI_REAL, MPI_SUM, comm2d, ierr )
377          IF ( cloud_physics ) THEN
378             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
379             CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb, &
380                                 MPI_REAL, MPI_SUM, comm2d, ierr )
381             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
382             CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb, &
383                                 MPI_REAL, MPI_SUM, comm2d, ierr )
384          ENDIF
385       ENDIF
386
387       IF ( passive_scalar )  THEN
388          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
389          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, &
390                              MPI_REAL, MPI_SUM, comm2d, ierr )
391       ENDIF
392#else
393       sums(:,1) = sums_l(:,1,0)
394       sums(:,2) = sums_l(:,2,0)
395       sums(:,4) = sums_l(:,4,0)
396       IF ( ocean )  sums(:,23) = sums_l(:,23,0)
397       IF ( humidity ) THEN
398          sums(:,44) = sums_l(:,44,0)
399          sums(:,41) = sums_l(:,41,0)
400          IF ( cloud_physics ) THEN
401             sums(:,42) = sums_l(:,42,0)
402             sums(:,43) = sums_l(:,43,0)
403          ENDIF
404       ENDIF
405       IF ( passive_scalar )  sums(:,41) = sums_l(:,41,0)
406#endif
407
408!
409!--    Final values are obtained by division by the total number of grid points
410!--    used for summation. After that store profiles.
411       sums(:,1) = sums(:,1) / ngp_2dh(sr)
412       sums(:,2) = sums(:,2) / ngp_2dh(sr)
413       sums(:,4) = sums(:,4) / ngp_2dh_s_inner(:,sr)
414       hom(:,1,1,sr) = sums(:,1)             ! u
415       hom(:,1,2,sr) = sums(:,2)             ! v
416       hom(:,1,4,sr) = sums(:,4)             ! pt
417
418
419!
420!--    Salinity
421       IF ( ocean )  THEN
422          sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr)
423          hom(:,1,23,sr) = sums(:,23)             ! sa
424       ENDIF
425
426!
427!--    Humidity and cloud parameters
428       IF ( humidity ) THEN
429          sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr)
430          sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
431          hom(:,1,44,sr) = sums(:,44)             ! vpt
432          hom(:,1,41,sr) = sums(:,41)             ! qv (q)
433          IF ( cloud_physics ) THEN
434             sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr)
435             sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr)
436             hom(:,1,42,sr) = sums(:,42)             ! qv
437             hom(:,1,43,sr) = sums(:,43)             ! pt
438          ENDIF
439       ENDIF
440
441!
442!--    Passive scalar
443       IF ( passive_scalar )  hom(:,1,41,sr) = sums(:,41) /  &
444            ngp_2dh_s_inner(:,sr)                    ! s (q)
445
446!
447!--    Horizontally averaged profiles of the remaining prognostic variables,
448!--    variances, the total and the perturbation energy (single values in last
449!--    column of sums_l) and some diagnostic quantities.
450!--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
451!--    ----  speaking the following k-loop would have to be split up and
452!--          rearranged according to the staggered grid.
453!--          However, this implies no error since staggered velocity components
454!--          are zero at the walls and inside buildings.
455       tn = 0
456#if defined( __intel_openmp_bug )
457       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, &
458       !$OMP                    tn, ust, ust2, u2, vst, vst2, v2, w2 )
459       tn = omp_get_thread_num()
460#else
461       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2, w2 )
462!$     tn = omp_get_thread_num()
463#endif
464       !$OMP DO
465       DO  i = nxl, nxr
466          DO  j =  nys, nyn
467             sums_l_etot = 0.0
468             DO  k = nzb_s_inner(j,i), nzt+1
469!
470!--             Prognostic and diagnostic variables
471                sums_l(k,3,tn)  = sums_l(k,3,tn)  + w(k,j,i)  * rmask(j,i,sr)
472                sums_l(k,8,tn)  = sums_l(k,8,tn)  + e(k,j,i)  * rmask(j,i,sr)
473                sums_l(k,9,tn)  = sums_l(k,9,tn)  + km(k,j,i) * rmask(j,i,sr)
474                sums_l(k,10,tn) = sums_l(k,10,tn) + kh(k,j,i) * rmask(j,i,sr)
475                sums_l(k,40,tn) = sums_l(k,40,tn) + p(k,j,i)
476
477                sums_l(k,33,tn) = sums_l(k,33,tn) + &
478                                  ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr)
479
480                IF ( humidity )  THEN
481                   sums_l(k,70,tn) = sums_l(k,70,tn) + &
482                                  ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr)
483                ENDIF
484
485!
486!--             Higher moments
487!--             (Computation of the skewness of w further below)
488                sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i)**3 * rmask(j,i,sr)
489
490                sums_l_etot  = sums_l_etot + &
491                                        0.5 * ( u(k,j,i)**2 + v(k,j,i)**2 +    &
492                                        w(k,j,i)**2 ) * rmask(j,i,sr)
493             ENDDO
494!
495!--          Total and perturbation energy for the total domain (being
496!--          collected in the last column of sums_l). Summation of these
497!--          quantities is seperated from the previous loop in order to
498!--          allow vectorization of that loop.
499             sums_l(nzb+4,pr_palm,tn) = sums_l(nzb+4,pr_palm,tn) + sums_l_etot
500!
501!--          2D-arrays (being collected in the last column of sums_l)
502             sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +   &
503                                        us(j,i)   * rmask(j,i,sr)
504             sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) + &
505                                        usws(j,i) * rmask(j,i,sr)
506             sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) + &
507                                        vsws(j,i) * rmask(j,i,sr)
508             sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) + &
509                                        ts(j,i)   * rmask(j,i,sr)
510             IF ( humidity )  THEN
511                sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) + &
512                                            qs(j,i)   * rmask(j,i,sr)
513             ENDIF
514          ENDDO
515       ENDDO
516
517!
518!--    Computation of statistics when ws-scheme is not used. Else these
519!--    quantities are evaluated in the advection routines.
520       IF ( .NOT. ws_scheme_mom .OR. sr /= 0 )  THEN
521          !$OMP DO
522          DO  i = nxl, nxr
523             DO  j =  nys, nyn
524                sums_l_eper = 0.0
525                DO  k = nzb_s_inner(j,i), nzt+1
526                   u2   = u(k,j,i)**2
527                   v2   = v(k,j,i)**2
528                   w2   = w(k,j,i)**2
529                   ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
530                   vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
531
532                   sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr)
533                   sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr)
534                   sums_l(k,32,tn) = sums_l(k,32,tn) + w2   * rmask(j,i,sr)
535!
536!--             Perturbation energy
537
538                   sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5 *       &
539                                  ( ust2 + vst2 + w2 ) * rmask(j,i,sr)
540                   sums_l_eper  = sums_l_eper + &
541                                        0.5 * ( ust2+vst2+w2 ) * rmask(j,i,sr)
542
543                ENDDO
544                sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn)   &
545                     + sums_l_eper
546             ENDDO
547          ENDDO
548       ENDIF
549!
550!--    Horizontally averaged profiles of the vertical fluxes
551
552       !$OMP DO
553       DO  i = nxl, nxr
554          DO  j = nys, nyn
555!
556!--          Subgridscale fluxes (without Prandtl layer from k=nzb,
557!--          oterwise from k=nzb+1)
558!--          NOTE: for simplicity, nzb_diff_s_inner is used below, although
559!--          ----  strictly speaking the following k-loop would have to be
560!--                split up according to the staggered grid.
561!--                However, this implies no error since staggered velocity
562!--                components are zero at the walls and inside buildings.
563
564             DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
565!
566!--             Momentum flux w"u"
567                sums_l(k,12,tn) = sums_l(k,12,tn) - 0.25 * (                   &
568                               km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
569                                                           ) * (               &
570                                   ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
571                                 + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
572                                                               ) * rmask(j,i,sr)
573!
574!--             Momentum flux w"v"
575                sums_l(k,14,tn) = sums_l(k,14,tn) - 0.25 * (                   &
576                               km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
577                                                           ) * (               &
578                                   ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
579                                 + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
580                                                               ) * rmask(j,i,sr)
581!
582!--             Heat flux w"pt"
583                sums_l(k,16,tn) = sums_l(k,16,tn)                              &
584                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
585                                               * ( pt(k+1,j,i) - pt(k,j,i) )   &
586                                               * ddzu(k+1) * rmask(j,i,sr)
587
588
589!
590!--             Salinity flux w"sa"
591                IF ( ocean )  THEN
592                   sums_l(k,65,tn) = sums_l(k,65,tn)                           &
593                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
594                                               * ( sa(k+1,j,i) - sa(k,j,i) )   &
595                                               * ddzu(k+1) * rmask(j,i,sr)
596                ENDIF
597
598!
599!--             Buoyancy flux, water flux (humidity flux) w"q"
600                IF ( humidity ) THEN
601                   sums_l(k,45,tn) = sums_l(k,45,tn)                           &
602                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
603                                               * ( vpt(k+1,j,i) - vpt(k,j,i) ) &
604                                               * ddzu(k+1) * rmask(j,i,sr)
605                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
606                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
607                                               * ( q(k+1,j,i) - q(k,j,i) )     &
608                                               * ddzu(k+1) * rmask(j,i,sr)
609
610                   IF ( cloud_physics ) THEN
611                      sums_l(k,51,tn) = sums_l(k,51,tn)                        &
612                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
613                                               * ( ( q(k+1,j,i) - ql(k+1,j,i) )&
614                                                - ( q(k,j,i) - ql(k,j,i) ) )   &
615                                               * ddzu(k+1) * rmask(j,i,sr) 
616                   ENDIF
617                ENDIF
618
619!
620!--             Passive scalar flux
621                IF ( passive_scalar )  THEN
622                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
623                                         - 0.5 * ( kh(k,j,i) + kh(k+1,j,i) )   &
624                                               * ( q(k+1,j,i) - q(k,j,i) )     &
625                                               * ddzu(k+1) * rmask(j,i,sr)
626                ENDIF
627
628             ENDDO
629
630!
631!--          Subgridscale fluxes in the Prandtl layer
632             IF ( use_surface_fluxes )  THEN
633                sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
634                                    usws(j,i) * rmask(j,i,sr)     ! w"u"
635                sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
636                                    vsws(j,i) * rmask(j,i,sr)     ! w"v"
637                sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
638                                    shf(j,i)  * rmask(j,i,sr)     ! w"pt"
639                sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
640                                    0.0 * rmask(j,i,sr)           ! u"pt"
641                sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
642                                    0.0 * rmask(j,i,sr)           ! v"pt"
643                IF ( ocean )  THEN
644                   sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
645                                       saswsb(j,i) * rmask(j,i,sr)  ! w"sa"
646                ENDIF
647                IF ( humidity )  THEN
648                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
649                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
650                   sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (           &
651                                       ( 1.0 + 0.61 * q(nzb,j,i) ) *   &
652                                       shf(j,i) + 0.61 * pt(nzb,j,i) * &
653                                                  qsws(j,i) )
654                   IF ( cloud_droplets )  THEN
655                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (        &
656                                         ( 1.0 + 0.61 * q(nzb,j,i) -   &
657                                           ql(nzb,j,i) ) * shf(j,i) +  &
658                                           0.61 * pt(nzb,j,i) * qsws(j,i) )
659                   ENDIF
660                   IF ( cloud_physics )  THEN
661!
662!--                   Formula does not work if ql(nzb) /= 0.0
663                      sums_l(nzb,51,tn) = sums_l(nzb,51,tn) + &   ! w"q" (w"qv")
664                                          qsws(j,i) * rmask(j,i,sr)
665                   ENDIF
666                ENDIF
667                IF ( passive_scalar )  THEN
668                   sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + &
669                                       qsws(j,i) * rmask(j,i,sr)  ! w"q" (w"qv")
670                ENDIF
671             ENDIF
672
673!
674!--          Subgridscale fluxes at the top surface
675             IF ( use_top_fluxes )  THEN
676                sums_l(nzt:nzt+1,12,tn) = sums_l(nzt:nzt+1,12,tn) + &
677                                    uswst(j,i) * rmask(j,i,sr)    ! w"u"
678                sums_l(nzt:nzt+1,14,tn) = sums_l(nzt:nzt+1,14,tn) + &
679                                    vswst(j,i) * rmask(j,i,sr)    ! w"v"
680                sums_l(nzt:nzt+1,16,tn) = sums_l(nzt:nzt+1,16,tn) + &
681                                    tswst(j,i)  * rmask(j,i,sr)   ! w"pt"
682                sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + &
683                                    0.0 * rmask(j,i,sr)           ! u"pt"
684                sums_l(nzt:nzt+1,61,tn) = sums_l(nzt:nzt+1,61,tn) + &
685                                    0.0 * rmask(j,i,sr)           ! v"pt"
686
687                IF ( ocean )  THEN
688                   sums_l(nzt,65,tn) = sums_l(nzt,65,tn) + &
689                                       saswst(j,i) * rmask(j,i,sr)  ! w"sa"
690                ENDIF
691                IF ( humidity )  THEN
692                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
693                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
694                   sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (             &
695                                       ( 1.0 + 0.61 * q(nzt,j,i) ) *     &
696                                       tswst(j,i) + 0.61 * pt(nzt,j,i) * &
697                                                  qswst(j,i) )
698                   IF ( cloud_droplets )  THEN
699                      sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (           &
700                                          ( 1.0 + 0.61 * q(nzt,j,i) -     &
701                                            ql(nzt,j,i) ) * tswst(j,i) +  &
702                                            0.61 * pt(nzt,j,i) * qswst(j,i) )
703                   ENDIF
704                   IF ( cloud_physics )  THEN
705!
706!--                   Formula does not work if ql(nzb) /= 0.0
707                      sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + &   ! w"q" (w"qv")
708                                          qswst(j,i) * rmask(j,i,sr)
709                   ENDIF
710                ENDIF
711                IF ( passive_scalar )  THEN
712                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + &
713                                       qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv")
714                ENDIF
715             ENDIF
716
717!
718!--          Resolved fluxes (can be computed for all horizontal points)
719!--          NOTE: for simplicity, nzb_s_inner is used below, although strictly
720!--          ----  speaking the following k-loop would have to be split up and
721!--                rearranged according to the staggered grid.
722             DO  k = nzb_s_inner(j,i), nzt
723                ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
724                              u(k+1,j,i) - hom(k+1,1,1,sr) )
725                vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
726                              v(k+1,j,i) - hom(k+1,1,2,sr) )
727                pts = 0.5 * ( pt(k,j,i)   - hom(k,1,4,sr) + &
728                              pt(k+1,j,i) - hom(k+1,1,4,sr) )
729
730!--             Higher moments
731                sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 * &
732                                                    rmask(j,i,sr)
733                sums_l(k,36,tn) = sums_l(k,36,tn) + pts**2 * w(k,j,i) * &
734                                                    rmask(j,i,sr)
735
736!
737!--             Salinity flux and density (density does not belong to here,
738!--             but so far there is no other suitable place to calculate)
739                IF ( ocean )  THEN
740                   IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
741                      pts = 0.5 * ( sa(k,j,i)   - hom(k,1,23,sr) + &
742                                 sa(k+1,j,i) - hom(k+1,1,23,sr) )
743                      sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) * &
744                                                       rmask(j,i,sr)
745                   ENDIF
746                   sums_l(k,64,tn) = sums_l(k,64,tn) + rho(k,j,i) * &
747                                                       rmask(j,i,sr)
748                   sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) * &
749                                                       rmask(j,i,sr)
750                ENDIF
751
752!
753!--             Buoyancy flux, water flux, humidity flux, liquid water
754!--             content, rain drop concentration and rain water content
755                IF ( humidity )  THEN
756                   IF ( cloud_physics .OR. cloud_droplets )  THEN
757                      pts = 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) +        &
758                                    vpt(k+1,j,i) - hom(k+1,1,44,sr) )
759                      sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
760                                                          rmask(j,i,sr)
761                      IF ( .NOT. cloud_droplets )  THEN
762                         pts = 0.5 *                                          &
763                              ( ( q(k,j,i) - ql(k,j,i) ) -                    &
764                              hom(k,1,42,sr) +                                &
765                              ( q(k+1,j,i) - ql(k+1,j,i) ) -                  &
766                              hom(k+1,1,42,sr) )
767                         sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) * &
768                                                             rmask(j,i,sr)
769                         IF ( icloud_scheme == 0  )  THEN
770                            sums_l(k,54,tn) = sums_l(k,54,tn) + ( ql(k,j,i) + &
771                                                                qr(k,j,i) ) * &
772                                                                rmask(j,i,sr)
773                            sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) *   &
774                                                                rmask(j,i,sr)
775                            sums_l(k,74,tn) = sums_l(k,74,tn) + qr(k,j,i) *   &
776                                                                rmask(j,i,sr)
777                            sums_l(k,75,tn) = sums_l(k,75,tn) + ql(k,j,i) *   &
778                                                                rmask(j,i,sr)
779                            sums_l(k,76,tn) = sums_l(k,76,tn) + prr(k,j,i) *  &
780                                                                rmask(j,i,sr)
781                         ELSE
782                            sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) *   &
783                                                                rmask(j,i,sr)
784                         ENDIF
785                      ELSE
786                         sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) *      &
787                                                             rmask(j,i,sr)
788                      ENDIF
789                   ELSE
790                      IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
791                         pts = 0.5 * ( vpt(k,j,i)   - hom(k,1,44,sr) + &
792                                       vpt(k+1,j,i) - hom(k+1,1,44,sr) )
793                         sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * &
794                                                             rmask(j,i,sr)
795                      ELSE IF ( ws_scheme_sca .AND. sr == 0 )  THEN
796                         sums_l(k,46,tn) = ( 1.0 + 0.61 * hom(k,1,41,sr) ) *  &
797                                             sums_l(k,17,tn) +      &
798                                          0.61 * hom(k,1,4,sr) * sums_l(k,49,tn)
799                      END IF
800                   END IF
801                ENDIF
802!
803!--             Passive scalar flux
804                IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca   &
805                     .OR. sr /= 0 ) )  THEN
806                   pts = 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) + &
807                                 q(k+1,j,i) - hom(k+1,1,41,sr) )
808                   sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
809                                                       rmask(j,i,sr)
810                ENDIF
811
812!
813!--             Energy flux w*e*
814!--             has to be adjusted
815                sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5 *          &
816                                             ( ust**2 + vst**2 + w(k,j,i)**2 )&
817                                             * rmask(j,i,sr)
818             ENDDO
819          ENDDO
820       ENDDO
821!
822!--    For speed optimization fluxes which have been computed in part directly
823!--    inside the WS advection routines are treated seperatly
824!--    Momentum fluxes first:
825       IF ( .NOT. ws_scheme_mom .OR. sr /= 0  )  THEN
826         !$OMP DO
827         DO  i = nxl, nxr
828            DO  j = nys, nyn
829               DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
830                  ust = 0.5 * ( u(k,j,i)   - hom(k,1,1,sr) + &
831                              u(k+1,j,i) - hom(k+1,1,1,sr) )
832                  vst = 0.5 * ( v(k,j,i)   - hom(k,1,2,sr) + &
833                              v(k+1,j,i) - hom(k+1,1,2,sr) )
834!
835!--               Momentum flux w*u*
836                  sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5 *                   &
837                                                    ( w(k,j,i-1) + w(k,j,i) ) &
838                                                    * ust * rmask(j,i,sr)
839!
840!--               Momentum flux w*v*
841                  sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5 *                   &
842                                                    ( w(k,j-1,i) + w(k,j,i) ) &
843                                                    * vst * rmask(j,i,sr)
844               ENDDO
845            ENDDO
846         ENDDO
847
848       ENDIF
849       IF ( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
850         !$OMP DO
851         DO  i = nxl, nxr
852            DO  j = nys, nyn
853               DO  k = nzb_diff_s_inner(j,i)-1, nzt_diff
854!
855!--               Vertical heat flux
856                  sums_l(k,17,tn) = sums_l(k,17,tn) +  0.5 * &
857                           ( pt(k,j,i)   - hom(k,1,4,sr) + &
858                           pt(k+1,j,i) - hom(k+1,1,4,sr) ) &
859                           * w(k,j,i) * rmask(j,i,sr)
860                  IF ( humidity )  THEN
861                     pts = 0.5 * ( q(k,j,i)   - hom(k,1,41,sr) + &
862                                q(k+1,j,i) - hom(k+1,1,41,sr) )
863                     sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * &
864                                                      rmask(j,i,sr)
865                  ENDIF
866               ENDDO
867            ENDDO
868         ENDDO
869
870       ENDIF
871
872
873!
874!--    Density at top follows Neumann condition
875       IF ( ocean )  THEN
876          sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
877          sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
878       ENDIF
879
880!
881!--    Divergence of vertical flux of resolved scale energy and pressure
882!--    fluctuations as well as flux of pressure fluctuation itself (68).
883!--    First calculate the products, then the divergence.
884!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
885       IF ( hom(nzb+1,2,55,0) /= 0.0  .OR.  hom(nzb+1,2,68,0) /= 0.0 )  THEN
886
887          sums_ll = 0.0  ! local array
888
889          !$OMP DO
890          DO  i = nxl, nxr
891             DO  j = nys, nyn
892                DO  k = nzb_s_inner(j,i)+1, nzt
893
894                   sums_ll(k,1) = sums_ll(k,1) + 0.5 * w(k,j,i) * (        &
895                  ( 0.25 * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1)   &
896                           - 0.5 * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) )     &
897                           ) )**2                                          &
898                + ( 0.25 * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i)   &
899                           - 0.5 * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) )     &
900                           ) )**2                                          &
901                   + w(k,j,i)**2                                  )
902
903                   sums_ll(k,2) = sums_ll(k,2) + 0.5 * w(k,j,i) &
904                                               * ( p(k,j,i) + p(k+1,j,i) )
905
906                ENDDO
907             ENDDO
908          ENDDO
909          sums_ll(0,1)     = 0.0    ! because w is zero at the bottom
910          sums_ll(nzt+1,1) = 0.0
911          sums_ll(0,2)     = 0.0
912          sums_ll(nzt+1,2) = 0.0
913
914          DO  k = nzb+1, nzt
915             sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
916             sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
917             sums_l(k,68,tn) = sums_ll(k,2)
918          ENDDO
919          sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
920          sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
921          sums_l(nzb,68,tn) = 0.0    ! because w* = 0 at nzb
922
923       ENDIF
924
925!
926!--    Divergence of vertical flux of SGS TKE and the flux itself (69)
927       IF ( hom(nzb+1,2,57,0) /= 0.0  .OR.  hom(nzb+1,2,69,0) /= 0.0 )  THEN
928
929          !$OMP DO
930          DO  i = nxl, nxr
931             DO  j = nys, nyn
932                DO  k = nzb_s_inner(j,i)+1, nzt
933
934                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5 * (                 &
935                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
936                 - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
937                                                             ) * ddzw(k)
938
939                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5 * (                 &
940                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
941                                                              )
942
943                ENDDO
944             ENDDO
945          ENDDO
946          sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
947          sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
948
949       ENDIF
950
951!
952!--    Horizontal heat fluxes (subgrid, resolved, total).
953!--    Do it only, if profiles shall be plotted.
954       IF ( hom(nzb+1,2,58,0) /= 0.0 ) THEN
955
956          !$OMP DO
957          DO  i = nxl, nxr
958             DO  j = nys, nyn
959                DO  k = nzb_s_inner(j,i)+1, nzt
960!
961!--                Subgrid horizontal heat fluxes u"pt", v"pt"
962                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5 *                   &
963                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
964                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
965                                                 * ddx * rmask(j,i,sr)
966                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5 *                   &
967                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
968                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
969                                                 * ddy * rmask(j,i,sr)
970!
971!--                Resolved horizontal heat fluxes u*pt*, v*pt*
972                   sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
973                                                  ( u(k,j,i) - hom(k,1,1,sr) ) &
974                                       * 0.5 * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
975                                                 pt(k,j,i)   - hom(k,1,4,sr) )
976                   pts = 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
977                                 pt(k,j,i)   - hom(k,1,4,sr) )
978                   sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
979                                                  ( v(k,j,i) - hom(k,1,2,sr) ) &
980                                       * 0.5 * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
981                                                 pt(k,j,i)   - hom(k,1,4,sr) )
982                ENDDO
983             ENDDO
984          ENDDO
985!
986!--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
987          sums_l(nzb,58,tn) = 0.0
988          sums_l(nzb,59,tn) = 0.0
989          sums_l(nzb,60,tn) = 0.0
990          sums_l(nzb,61,tn) = 0.0
991          sums_l(nzb,62,tn) = 0.0
992          sums_l(nzb,63,tn) = 0.0
993
994       ENDIF
995
996!
997!--    Calculate the user-defined profiles
998       CALL user_statistics( 'profiles', sr, tn )
999       !$OMP END PARALLEL
1000
1001!
1002!--    Summation of thread sums
1003       IF ( threads_per_task > 1 )  THEN
1004          DO  i = 1, threads_per_task-1
1005             sums_l(:,3,0)          = sums_l(:,3,0) + sums_l(:,3,i)
1006             sums_l(:,4:40,0)       = sums_l(:,4:40,0) + sums_l(:,4:40,i)
1007             sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
1008                                      sums_l(:,45:pr_palm,i)
1009             IF ( max_pr_user > 0 )  THEN
1010                sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
1011                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
1012                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
1013             ENDIF
1014          ENDDO
1015       ENDIF
1016
1017#if defined( __parallel )
1018
1019!
1020!--    Compute total sum from local sums
1021       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1022       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, &
1023                           MPI_SUM, comm2d, ierr )
1024#else
1025       sums = sums_l(:,:,0)
1026#endif
1027
1028!
1029!--    Final values are obtained by division by the total number of grid points
1030!--    used for summation. After that store profiles.
1031!--    Profiles:
1032       DO  k = nzb, nzt+1
1033          sums(k,3)               = sums(k,3)           / ngp_2dh(sr)
1034          sums(k,8:11)            = sums(k,8:11)        / ngp_2dh_s_inner(k,sr)
1035          sums(k,12:22)           = sums(k,12:22)       / ngp_2dh(sr)
1036          sums(k,23:29)           = sums(k,23:29)       / ngp_2dh_s_inner(k,sr)
1037          sums(k,30:32)           = sums(k,30:32)       / ngp_2dh(sr)
1038          sums(k,33:34)           = sums(k,33:34)       / ngp_2dh_s_inner(k,sr)
1039          sums(k,35:39)           = sums(k,35:39)       / ngp_2dh(sr)
1040          sums(k,40)              = sums(k,40)          / ngp_2dh_s_inner(k,sr)
1041          sums(k,45:53)           = sums(k,45:53)       / ngp_2dh(sr)
1042          sums(k,54)              = sums(k,54)          / ngp_2dh_s_inner(k,sr)
1043          sums(k,55:63)           = sums(k,55:63)       / ngp_2dh(sr)
1044          sums(k,64)              = sums(k,64)          / ngp_2dh_s_inner(k,sr)
1045          sums(k,65:69)           = sums(k,65:69)       / ngp_2dh(sr)
1046          sums(k,70:pr_palm-2)    = sums(k,70:pr_palm-2)/ ngp_2dh_s_inner(k,sr)
1047       ENDDO
1048
1049!--    Upstream-parts
1050       sums(nzb:nzb+11,pr_palm-1) = sums(nzb:nzb+11,pr_palm-1) / ngp_3d(sr)
1051!--    u* and so on
1052!--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
1053!--    size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
1054!--    above the topography, they are being divided by ngp_2dh(sr)
1055       sums(nzb:nzb+3,pr_palm)    = sums(nzb:nzb+3,pr_palm)    / &
1056                                    ngp_2dh(sr)
1057       sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)       / &    ! qs
1058                                    ngp_2dh(sr)
1059!--    eges, e*
1060       sums(nzb+4:nzb+5,pr_palm)  = sums(nzb+4:nzb+5,pr_palm)  / &
1061                                    ngp_3d(sr)
1062!--    Old and new divergence
1063       sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
1064                                    ngp_3d_inner(sr)
1065
1066!--    User-defined profiles
1067       IF ( max_pr_user > 0 )  THEN
1068          DO  k = nzb, nzt+1
1069             sums(k,pr_palm+1:pr_palm+max_pr_user) = &
1070                                    sums(k,pr_palm+1:pr_palm+max_pr_user) / &
1071                                    ngp_2dh_s_inner(k,sr)
1072          ENDDO
1073       ENDIF
1074
1075!
1076!--    Collect horizontal average in hom.
1077!--    Compute deduced averages (e.g. total heat flux)
1078       hom(:,1,3,sr)  = sums(:,3)      ! w
1079       hom(:,1,8,sr)  = sums(:,8)      ! e     profiles 5-7 are initial profiles
1080       hom(:,1,9,sr)  = sums(:,9)      ! km
1081       hom(:,1,10,sr) = sums(:,10)     ! kh
1082       hom(:,1,11,sr) = sums(:,11)     ! l
1083       hom(:,1,12,sr) = sums(:,12)     ! w"u"
1084       hom(:,1,13,sr) = sums(:,13)     ! w*u*
1085       hom(:,1,14,sr) = sums(:,14)     ! w"v"
1086       hom(:,1,15,sr) = sums(:,15)     ! w*v*
1087       hom(:,1,16,sr) = sums(:,16)     ! w"pt"
1088       hom(:,1,17,sr) = sums(:,17)     ! w*pt*
1089       hom(:,1,18,sr) = sums(:,16) + sums(:,17)    ! wpt
1090       hom(:,1,19,sr) = sums(:,12) + sums(:,13)    ! wu
1091       hom(:,1,20,sr) = sums(:,14) + sums(:,15)    ! wv
1092       hom(:,1,21,sr) = sums(:,21)     ! w*pt*BC
1093       hom(:,1,22,sr) = sums(:,16) + sums(:,21)    ! wptBC
1094                                       ! profile 24 is initial profile (sa)
1095                                       ! profiles 25-29 left empty for initial
1096                                       ! profiles
1097       hom(:,1,30,sr) = sums(:,30)     ! u*2
1098       hom(:,1,31,sr) = sums(:,31)     ! v*2
1099       hom(:,1,32,sr) = sums(:,32)     ! w*2
1100       hom(:,1,33,sr) = sums(:,33)     ! pt*2
1101       hom(:,1,34,sr) = sums(:,34)     ! e*
1102       hom(:,1,35,sr) = sums(:,35)     ! w*2pt*
1103       hom(:,1,36,sr) = sums(:,36)     ! w*pt*2
1104       hom(:,1,37,sr) = sums(:,37)     ! w*e*
1105       hom(:,1,38,sr) = sums(:,38)     ! w*3
1106       hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20 )**1.5   ! Sw
1107       hom(:,1,40,sr) = sums(:,40)     ! p
1108       hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
1109       hom(:,1,46,sr) = sums(:,46)     ! w*vpt*       
1110       hom(:,1,47,sr) = sums(:,45) + sums(:,46)    ! wvpt
1111       hom(:,1,48,sr) = sums(:,48)     ! w"q" (w"qv")
1112       hom(:,1,49,sr) = sums(:,49)     ! w*q* (w*qv*)
1113       hom(:,1,50,sr) = sums(:,48) + sums(:,49)    ! wq (wqv)
1114       hom(:,1,51,sr) = sums(:,51)     ! w"qv"
1115       hom(:,1,52,sr) = sums(:,52)     ! w*qv*       
1116       hom(:,1,53,sr) = sums(:,52) + sums(:,51)    ! wq (wqv)
1117       hom(:,1,54,sr) = sums(:,54)     ! ql
1118       hom(:,1,55,sr) = sums(:,55)     ! w*u*u*/dz
1119       hom(:,1,56,sr) = sums(:,56)     ! w*p*/dz
1120       hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho )/dz
1121       hom(:,1,58,sr) = sums(:,58)     ! u"pt"
1122       hom(:,1,59,sr) = sums(:,59)     ! u*pt*
1123       hom(:,1,60,sr) = sums(:,58) + sums(:,59)    ! upt_t
1124       hom(:,1,61,sr) = sums(:,61)     ! v"pt"
1125       hom(:,1,62,sr) = sums(:,62)     ! v*pt*
1126       hom(:,1,63,sr) = sums(:,61) + sums(:,62)    ! vpt_t
1127       hom(:,1,64,sr) = sums(:,64)     ! rho
1128       hom(:,1,65,sr) = sums(:,65)     ! w"sa"
1129       hom(:,1,66,sr) = sums(:,66)     ! w*sa*
1130       hom(:,1,67,sr) = sums(:,65) + sums(:,66)    ! wsa
1131       hom(:,1,68,sr) = sums(:,68)     ! w*p*
1132       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho
1133       hom(:,1,70,sr) = sums(:,70)     ! q*2
1134       hom(:,1,71,sr) = sums(:,71)     ! prho
1135       hom(:,1,72,sr) = hyp * 1E-4     ! hyp in dbar
1136       hom(:,1,73,sr) = sums(:,73)     ! nr
1137       hom(:,1,74,sr) = sums(:,74)     ! qr
1138       hom(:,1,75,sr) = sums(:,75)     ! qc
1139       hom(:,1,76,sr) = sums(:,76)     ! prr (precipitation rate)
1140
1141       hom(:,1,pr_palm-1,sr) = sums(:,pr_palm-1)
1142                                       ! upstream-parts u_x, u_y, u_z, v_x,
1143                                       ! v_y, usw. (in last but one profile)
1144       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
1145                                       ! u*, w'u', w'v', t* (in last profile)
1146
1147       IF ( max_pr_user > 0 )  THEN    ! user-defined profiles
1148          hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
1149                               sums(:,pr_palm+1:pr_palm+max_pr_user)
1150       ENDIF
1151
1152!
1153!--    Determine the boundary layer height using two different schemes.
1154!--    First scheme: Starting from the Earth's (Ocean's) surface, look for the
1155!--    first relative minimum (maximum) of the total heat flux.
1156!--    The corresponding height is assumed as the boundary layer height, if it
1157!--    is less than 1.5 times the height where the heat flux becomes negative
1158!--    (positive) for the first time.
1159       z_i(1) = 0.0
1160       first = .TRUE.
1161
1162       IF ( ocean )  THEN
1163          DO  k = nzt, nzb+1, -1
1164             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
1165                .AND. abs(hom(k,1,18,sr)) > 1.0E-8)  THEN
1166                first = .FALSE.
1167                height = zw(k)
1168             ENDIF
1169             IF ( hom(k,1,18,sr) < 0.0  .AND. &
1170                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
1171                  hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
1172                IF ( zw(k) < 1.5 * height )  THEN
1173                   z_i(1) = zw(k)
1174                ELSE
1175                   z_i(1) = height
1176                ENDIF
1177                EXIT
1178             ENDIF
1179          ENDDO
1180       ELSE
1181          DO  k = nzb, nzt-1
1182             IF ( first .AND. hom(k,1,18,sr) < 0.0 &
1183                .AND. abs(hom(k,1,18,sr)) > 1.0E-8 )  THEN
1184                first = .FALSE.
1185                height = zw(k)
1186             ENDIF
1187             IF ( hom(k,1,18,sr) < 0.0  .AND. & 
1188                  abs(hom(k,1,18,sr)) > 1.0E-8 .AND. &
1189                  hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
1190                IF ( zw(k) < 1.5 * height )  THEN
1191                   z_i(1) = zw(k)
1192                ELSE
1193                   z_i(1) = height
1194                ENDIF
1195                EXIT
1196             ENDIF
1197          ENDDO
1198       ENDIF
1199
1200!
1201!--    Second scheme: Gradient scheme from Sullivan et al. (1998), modified
1202!--    by Uhlenbrock(2006). The boundary layer height is the height with the
1203!--    maximal local temperature gradient: starting from the second (the last
1204!--    but one) vertical gridpoint, the local gradient must be at least
1205!--    0.2K/100m and greater than the next four gradients.
1206!--    WARNING: The threshold value of 0.2K/100m must be adjusted for the
1207!--             ocean case!
1208       z_i(2) = 0.0
1209       DO  k = nzb+1, nzt+1
1210          dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
1211       ENDDO
1212       dptdz_threshold = 0.2 / 100.0
1213
1214       IF ( ocean )  THEN
1215          DO  k = nzt+1, nzb+5, -1
1216             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
1217                  dptdz(k) > dptdz(k-1)  .AND.  dptdz(k) > dptdz(k-2)  .AND.  &
1218                  dptdz(k) > dptdz(k-3)  .AND.  dptdz(k) > dptdz(k-4) )  THEN
1219                z_i(2) = zw(k-1)
1220                EXIT
1221             ENDIF
1222          ENDDO
1223       ELSE
1224          DO  k = nzb+1, nzt-3
1225             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
1226                  dptdz(k) > dptdz(k+1)  .AND.  dptdz(k) > dptdz(k+2)  .AND.  &
1227                  dptdz(k) > dptdz(k+3)  .AND.  dptdz(k) > dptdz(k+4) )  THEN
1228                z_i(2) = zw(k-1)
1229                EXIT
1230             ENDIF
1231          ENDDO
1232       ENDIF
1233
1234       hom(nzb+6,1,pr_palm,sr) = z_i(1)
1235       hom(nzb+7,1,pr_palm,sr) = z_i(2)
1236
1237!
1238!--    Computation of both the characteristic vertical velocity and
1239!--    the characteristic convective boundary layer temperature.
1240!--    The horizontal average at nzb+1 is input for the average temperature.
1241       IF ( hom(nzb,1,18,sr) > 0.0 .AND. abs(hom(nzb,1,18,sr)) > 1.0E-8 &
1242           .AND.  z_i(1) /= 0.0 )  THEN
1243          hom(nzb+8,1,pr_palm,sr)  = ( g / hom(nzb+1,1,4,sr) * &
1244                                       hom(nzb,1,18,sr) *      &
1245                                       ABS( z_i(1) ) )**0.333333333
1246!--       so far this only works if Prandtl layer is used
1247          hom(nzb+11,1,pr_palm,sr) = hom(nzb,1,16,sr) / hom(nzb+8,1,pr_palm,sr)
1248       ELSE
1249          hom(nzb+8,1,pr_palm,sr)  = 0.0
1250          hom(nzb+11,1,pr_palm,sr) = 0.0
1251       ENDIF
1252
1253!
1254!--    Collect the time series quantities
1255       ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)     ! E
1256       ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)     ! E*
1257       ts_value(3,sr) = dt_3d
1258       ts_value(4,sr) = hom(nzb,1,pr_palm,sr)       ! u*
1259       ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr)     ! th*
1260       ts_value(6,sr) = u_max
1261       ts_value(7,sr) = v_max
1262       ts_value(8,sr) = w_max
1263       ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr)    ! new divergence
1264       ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr)    ! old Divergence
1265       ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr)    ! z_i(1)
1266       ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr)    ! z_i(2)
1267       ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr)    ! w*
1268       ts_value(14,sr) = hom(nzb,1,16,sr)           ! w'pt'   at k=0
1269       ts_value(15,sr) = hom(nzb+1,1,16,sr)         ! w'pt'   at k=1
1270       ts_value(16,sr) = hom(nzb+1,1,18,sr)         ! wpt     at k=1
1271       ts_value(17,sr) = hom(nzb,1,4,sr)            ! pt(0)
1272       ts_value(18,sr) = hom(nzb+1,1,4,sr)          ! pt(zp)
1273       ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr)    ! u'w'    at k=0
1274       ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)    ! v'w'    at k=0
1275       ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
1276
1277       IF ( ts_value(5,sr) /= 0.0 )  THEN
1278          ts_value(22,sr) = ts_value(4,sr)**2 / &
1279                            ( kappa * g * ts_value(5,sr) / ts_value(18,sr) ) ! L
1280       ELSE
1281          ts_value(22,sr) = 10000.0
1282       ENDIF
1283
1284       ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
1285!
1286!--    Calculate additional statistics provided by the user interface
1287       CALL user_statistics( 'time_series', sr, 0 )
1288
1289    ENDDO    ! loop of the subregions
1290
1291!
1292!-- If required, sum up horizontal averages for subsequent time averaging
1293    IF ( do_sum )  THEN
1294       IF ( average_count_pr == 0 )  hom_sum = 0.0
1295       hom_sum = hom_sum + hom(:,1,:,:)
1296       average_count_pr = average_count_pr + 1
1297       do_sum = .FALSE.
1298    ENDIF
1299
1300!
1301!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
1302!-- This flag is reset after each time step in time_integration.
1303    flow_statistics_called = .TRUE.
1304
1305    CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
1306
1307
1308 END SUBROUTINE flow_statistics
Note: See TracBrowser for help on using the repository browser.