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

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

New:
---

further adjustments on Tsubame and concerning openMP usage
(mrun, mbuild, subjob)

Changed:


Errors:


Bugfix in calculation of divergence of vertical flux of resolved scale
energy, pressure fluctuations, and flux of pressure fluctuation itself
(flow statistics)

Bugfix: module pegrid was missing. (user_statistics)

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