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

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