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

Last change on this file since 647 was 625, checked in by heinze, 13 years ago

last commit documented

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