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

Last change on this file since 132 was 132, checked in by letzel, 17 years ago

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