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

Last change on this file since 855 was 802, checked in by suehring, 12 years ago

Last commit documented.

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