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

Last change on this file since 133 was 133, checked in by letzel, 14 years ago

Vertical profiles now based on nzb_s_inner; they are divided by
ngp_2dh_s_inner (scalars, procucts of scalars) and ngp_2dh (staggered velocity
components and their products, procucts of scalars and velocity components),
respectively.

Allow new case bc_uv_t = 'dirichlet_0' for channel flow.

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