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

Last change on this file since 709 was 709, checked in by raasch, 10 years ago

formatting adjustments

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