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

Last change on this file since 1053 was 1053, checked in by hoffmann, 9 years ago

two-moment cloud physics implemented

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