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

Last change on this file since 1112 was 1112, checked in by raasch, 8 years ago

last commit documented

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