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

Last change on this file since 623 was 623, checked in by raasch, 14 years ago

last commit documented

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